      program disco
c	dust radiative transfer in a circumstellar disk
c	with nested r-theta grids tranformed onto rectangular grids
c	for the rad xfer
c	expand the innner grid by 8 to resolve the inner rim better
c	start with Chiang & Goldreich disk model
c	with hydro equil for vertical structure
c	then do full dust rad xfer

        implicit real*8 (a-h,o-z)
#include "mpif.h"
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        parameter (mc = 256)
        parameter (mmap = 1024)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,epspe,agas,sigma,
     &          colin,colup,xco,xwat
        common /angls/ dcost,dcos9,rot,irot
        common /parmz/ rgrid1,r0,sigma0,pcol,qcol,avisc,xmstar,xlstar,
     &          tuv,xuv,tsubl,gsubl,enhanc,tdiss,tice,xco0,xwat0,
     &          coswind,xc,xsi,ac,asi,kit,lit,mit,nnt,lx
        common /xj/ xju,xjp,xjq,xjr,xk
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)
        common /taus/ taucont,tauline,tauterm,iuv,iusil,iufl,iwterm
        common /cube/ cube(mc,mc,mc),ngc
        common /exptab/ expi(10239),too
        dimension mpistatus(MPI_STATUS_SIZE)
        dimension taud(mu,mt),dtac4(mt)
        dimension colr(mt,mr,mg),colz(mt,mr,mg),tauv(mt,mr,mg)
        dimension dt(mt,mr,mg),dm(mt,mr,mg)
        dimension fw(mw,ml,2*mt),fwj(mw)
        dimension spec(mw,mg)
        dimension frame(mc,mc),picture(mt,mr,mg)
        dimension con(mmap)
        logical ifmpi,doplot,doprint,wait
        logical smooth,lte
        logical ifrim
        logical badnum
        character*16 monfile,pfile,qfile,rfile,hfile,ifile,lockline

c	physical constants

        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0+00)
        parameter (ten = 1.0d+01)
        parameter (hun = 1.0d+02)
        parameter (thrq = 0.75d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (fifth = 0.2d+00)
        parameter (pi = 3.14159265d+00)
        parameter (sqrt2 = 1.414214d+00)
        parameter (tiny = 1.0d-12)
        parameter (zilch = 1.0d-200)
        parameter (clt = 2.9979d+10)
        parameter (hpl = 6.6262d-27)
        parameter (hc2 = 5.955d-06)
        parameter (bk = 1.3807d-16)
        parameter (sigsb = 5.670d-05)
        parameter (au = 1.496d+13)
        parameter (rsun = 6.96d+10)
        parameter (xlsun = 3.90d+33)
        parameter (xmsun = 1.99d+33)
c	mass / gas particle
        parameter (xmg = 4.0d-24)
        parameter (vkms = 1.0d+05)

        alog2 = log(two)
        sqrtpi = sqrt(pi)
        tpi = two*pi
        fpi = four*pi
        hc = hpl*clt
        hck = hc/bk

        ifmpi = .true.
        if(ifmpi) then
          call mpi_init(ierr)
          call mpi_comm_size(MPI_COMM_WORLD,mysize,ierr)
          call mpi_comm_rank(MPI_COMM_WORLD,myrank,ierr)
          ifmpi = (mysize.gt.1)
        else
          mysize = 1
          myrank = 0
        endif
        
        if(myrank.eq.mysize/2) then
          print '("Running disco with",i4," processors")', mysize
          doprint = .true.
          doplot = .true.
          print '(" Enter name of output file ",$)'
        else
          doprint = .false.
          doplot = .false.
        endif
        wait = .false.
        wait = .true.

c	initialize expi arrays
        x = expinit()

c	initialize plot files

        if(myrank.eq.0) then
          read '(a16)', monfile
          if(index(monfile,' ').lt.2) monfile = 'specs'
        endif
        if(ifmpi) call mpi_bcast(monfile,16,MPI_CHARACTER,
     &                           0,MPI_COMM_WORLD,ierr)
        if(doplot) then
          open(unit=8,file='plotlock',access='sequential')
          write(8,'("empty")')
          close(unit=8)
          last = index(monfile,' ')-1
          last = min(last,14)
          ifile = monfile(1:last)//'.log'
          open(unit=16,file=ifile,access='sequential')
          monfile = ''
        endif

c	input parameters

c	outer radius of grid 1
        rgrid1 = 16.*au
c	column density knee radius
        r0 = 10.0*au
c	gas column density in disk at r0 (g/cm2 converted to mol/cm2)
        sigma0 = 100.0/xmg
c	radial power law negative index for column density
c	C&G use 1.5, but D''A finds ~1.0 for constant alpha
        pcol = 1.0
c	radial power law positive index for column density inside r0
        qcol = 10.0
c	to specify sigma at 1 AU instead of inside r0
c        sigma0 = sigma0*((au/r0)**pcol+(r0/au)**qcol)
c	dust density and viscosity parameter for dust settling
        rhod = 3.0d+00
        avisc = 1.0d-02
c	radius of star
        rstar = 1.7*rsun
c	mass of star (M_sun)
        xmstar = 2.4
c	luminosity of star (L_sun)
        xlstar = 32.0
c	add a hot source of UV to stellar spectrum
        tuv = 1.5d+04
        xuv = 0.1
c	dust sublimation temperature coefficients (for both carbon and silcate)
c	tsubli = tsubl*pres**gsubl
        tsubl = 600.
        gsubl = 0.02
c	enhancement of stellar radiation by reflection and emission from disk
        enhanc = 0.6
c	molecule destruction temperature
c	CO should extend in to where T = 4000.
        tdiss = 1000.
c	molecule freeze-out temperature
        tice = 30.
c	accretion heating
        xmdot = 1.0d-08
c	dust to gas heat flow acommodation coefficient
        acom = 0.05d+00
c	x-ray luminosity of star
	pxr = zero
c	x-ray absorption cross section per H
c	~1e-22 at 1 keV, falling as ~E^-2.5 (Gorti&Hollenbach)
        axr = 1.0d-22
c	photoelectric heating of gas
        npe = nu-4
        epspe = 0.05d+00
c	abundances for cooling function; should vary in disk
        xco0 = 1.0d-04
        xwat0 = 1.0d-05
c	angular width of wind-swept cone
        coswind = 0.5
        thwind = acos(coswind)
        dvturb0 = 0.25*vkms
        dvturb1 = 0.25*vkms
c	carbon and silicon fraction in grains
        xc = 0.01
        xsi = 1.0
c	carbon and silicate rolloff radius (cm) (amax = 0.1 cm)
        ac = 1.0d-04
        asi = 1.0d-02
c	wavenumber (>0) or grain type (<0) for ggfs plots
        iupl = -6

c	C&G iterations
        kit = 4
c	dust radiative transfer iterations
        lit = 16
c	maximum line radiative transfer iterations
        mit = 16
c	number of grain sizes or types; first half are carbon
c	parameter (ma = 32)
        na = 24
c	number of wavelengths for dust spectrum
c	parameter (mu = 64)
        nu = 64
c	number of nested grids (<= 12)
c	parameter (mg = 20)
        ng = 12
c	grid density multiplier for inner grid
c	nn9 = 8 make pretty pictures, but it takes forever for the lines
        nn9 = 4
c	parameter (mr = 128)
        nr = 128
c	note: costh now runs from 1 to 0
c	set nnt = nr so r*dcosth = dr at outer edge of grids
c	or maybe bigger is better (and if mult of 3 can use 12 processors)
        nnt = 3*nr/2
c	but don''t bother calculating at theta above disk
c	there might be a problem unless nt/nn9 is a multiple of mysize
c	parameter (mt = 192)
        nt = nnt/2
c	nx = 16 seems to be OK for SED, but maybe not for disk structure
        nx = 64
c	use lx,lz,lrot for nx,nz,nrot for the line calculation
        lx = 16
c	make nrot a fraction of nnt (or vice versa?) and a multiple of mysize
        nrot = nnt/2
        if(mod(nrot,mysize).ne.0) then
          if(doprint) print '(" nrot,mysize =",2i6)', nrot,mysize
          stop
        endif
        lrot = nrot/2
        smooth = .true.
        if(smooth) then
c	angular smoothing of radiation field in grid 9.  ntsm = 2 in grids 1-8
          ntsm = 2
c	radial smoothing hardwired for nrsm = 2 if set > 1
          nrsm = 2
        else
          ntsm = 1
          nrsm = 1
        endif
	smooth = .false.  !  turn off extra angular smoothing in grid 9
c	outermost grid in cubes
        ngc = ng/2
c	grid for ggfs plots
        igpl = 4
c	parameter (mc = 256)
c	parameter (mmap = 1024)
        nc = mc

        call readparms

c	stellar effective temperature
        if(xlstar*rstar.gt.zero) then
          tstar = 5.8d+03*(xlstar/(rstar/rsun)**2)**quart
        elseif(xlstar*tstar.gt.zero) then
          rstar = rsun*sqrt(xlstar)*(5.8d+03/tstar)**2
        elseif(rstar*tstar.gt.zero) then
          xlstar = (tstar/5.8d+03)**4*(rstar/rsun)**2
        endif
c	T s.t. kT = G Mstar Mmol / Rstar
        tgrav = 6.67d-08*xmstar*xmsun*xmg/(rstar*bk)
c	effective temperature including UV excess
        teff = tstar*(one+xuv)**quart
c	hot dust temperature at rstar
        tsurf = teff/(two**0.4)
        gmmdot = 8.36d+51*xmstar*xmdot
c	viscosity parameter for dust settling
        avisc = avisc*xmg/rhod
c	x-ray luminosity of star * heating efficiency
        pxr = 0.4*pxr
c	orbital speed at 1 AU (scale as r^-1/2)
        vorb = 30.*vkms*sqrt(xmstar)
c	vterm may be hardwired = 0
c	and it may not work with bipolar hardwired
        vterm = zero*vkms
c	Gaussian resolution HW1/e in pixels
        dsm = vkms/dv
        nsm = 32
        ngr = ng+nn9
        xn9 = nn9
        xnn9 = xn9*xn9
        ng9 = ng+1
        nr9 = nn9*nr
c	cell size on grid 0; dr = dr0*2**ig
        dr0 = rgrid1/(2*nr)
c	cell size on expanded grid 1
        dr9 = two*dr0/xn9
        nr2 = nr/2
        nt8 = nt-nt/nn9
        dcost = one/nnt
        dcos9 = dcost/xn9
        cz8 = nt*dcos9
        coswind = min(coswind,(nt*dcost))
        thwind = acos(coswind)
        itwind = nt-int(coswind/dcost)
        dcrot = one/nrot
c	note: pisco now assumes bipolar symmetry
        if(vterm.eq.zero) then
          ly = lx/2
        else
          ly = lx
        endif
        nx9 = nn9*nx
        ny = nx
c	set nz = 2*nr so dz = dr
        nz = 2*max(nr,nnt)
        nz9 = nn9*nz
        xnx = nx
        xnz = nz

c	inner and outer edge of calculated disk
        rin = rgrid1/four
        rmin = r0*(tpi*(1.0d+03/tgrav)*r0**3/(sigma0**2*rstar))
     &            **(one/(two*qcol-3.))
        if((qcol.gt.two).and.(rin.lt.rmin)) rin = rmin
        rout = rgrid1*two**(ng-1)
        rmax = r0*(tpi*(1.0d+03/tgrav)*r0**3/(sigma0**2*rstar))
     &            **(one/(-two*pcol-3.))
        if((pcol.gt.zero).and.(rout.gt.rmax)) rout = rmax

c	pass parameters to startgr
        call initgr(dr0,dr9,dcost,dcos9,nr,nr9,ng9,nt,nt8)

        if((ngr.gt.mg).or.(nr.gt.mr).or.(nt.gt.mt).or.(nrot.gt.2*mt))
     &    then
          if(doprint) then
            print '(" array out of bounds",10i6,2i8)',
     &      ngr,mg,nr,mr,nt,mt,nrot,mt,nz,nt*ng*nr2,mmap
            stop
          endif
        endif
        if(ifmpi.and.doprint
     &    .and.(mod(nrot,mysize).ne.0)) then
          print '(" nrot must be a multiple of mysize",2i8)',
     &      nrot,mysize
          stop
        endif

        dx0 = (two*nr*dr0)/xnx
        dx9 = two*dx0/xn9
        dz0 = xnx*dx0/xnz
        dz9 = two*dz0/xn9
        nx2 = nx/2
        ny2 = ny/2
        nz2 = nz/2
        if(doprint) then
c	typical values
          tgi = 500.
          hh2 = (tgi/tgrav)*(r0**3/rstar)
          sigma = sigma0/two
          dens = sigma/sqrt(tpi*hh2)
          pm = two-pcol
          xmdisk = (tpi*r0**2/pm)*(sigma0*xmg/xmsun)*((rout/r0)**pm-one)
          print'(" disk mass, density, column density at",f8.3," AU =",
     &      1p,3e9.2)',
     &      r0/au,xmdisk,dens,sigma*xmg
          print '(" stellar M, L, T, R, dMdt =",
     &      2f8.3,f8.0,f8.2,1p,e9.2)',
     &      xmstar,xlstar,tstar,rstar/rsun,xmdot
          print '(" vorb, vturb, vstep =",f8.1,3f8.3)',
     &      vorb/vkms,dvturb0/vkms,dvturb1/vkms,dv/vkms
          print '(" inner xyz grid:",i3," x",f10.5," to",f10.4," AU")',
     &      nx9,dx9/au,nx9*dx9/(two*au)
          print '(" outer xyz grid:",i3," x",f10.4," to",f10.3," AU")',
     &      nx,(two**ng)*dx0/au,(two**(ng-1))*nx*dx0/au
          write(16,'(" disk mass, density, column density at",f8.3,
     &      " AU =",1p,3e9.2)')
     &      r0/au,xmdisk,dens,sigma*xmg
          write(16,'(" stellar M, L, T, R, dMdt =",
     &      2f8.3,f8.0,f8.2,1p,e9.2)')
     &      xmstar,xlstar,tstar,rstar/rsun,xmdot
          write(16,'(" vorb, vturb, vstep =",f8.1,3f8.3)')
     &      vorb/vkms,dvturb0/vkms,dvturb1/vkms,dv/vkms
          write(16,'(" inner xyz grid:",i3," x",f10.5,
     &      " to",f10.4," AU")')
     &       nx9,dx9/au,nx9*dx9/(two*au)
          write(16,'(" outer xyz grid:",i3," x",f10.4,
     &      " to",f10.3," AU")')
     &      nx,(two**ng)*dx0/au,(two**(ng-1))*nx*dx0/au
        endif
        astar = pi*(rstar/(dx0*two**ng))**2
        dxstar = rstar/four
        gsum = zero
        do ix = 1,4
          x = dxstar*(ix-half)
          gwt = sqrt(one-(x/rstar)**2)
          gsum = gsum+two*gwt
        enddo
        bstar = astar/gsum

c	set up dust arrays
        if(doprint) print '(" initializing dust property arrays")'
        tmid = 1.0d+03
        hh2 = (tmid/tgrav)*(r0**3/rstar)
        pres = tmid*sigma0/sqrt(tpi*hh2)
        tsubli = tsubl*max(one,pres)**gsubl
        call dustinit(xc,xsi,ac,asi,tuv,xuv,emiss,tsubli)
c	abss is the fraction of the incident starlight absorbed by the disk
        abss = sqrt(emiss)
        fwmx = zero
        sumfs = zero
        do iu = 1,nu
          busdud = bus(iu)*dud(iu)
          if(busdud.gt.fwmx) fwmx = busdud
          sumfs = sumfs+busdud
          gg(iu,1,1,1) = zero
          opa(iu,1,1,1) = zero
        enddo
        fwmx = astar*fwmx
        if(doprint) then
          xlvis = sigsb*tstar**4*fpi*rstar**2/xlsun
          xluv = xuv*xlvis
          xlphot = sumfs*pi*fpi*rstar**2/xlsun
          xlaccr = gmmdot/(two*rin*xlsun)
          print '(" Lvis,Luv,Lphot,Laccr(rin,rstar) =",5f8.3)',
     &      xlvis,xluv,xlphot,xlaccr,xlaccr*(rin/rstar)
c	heating and cooling power per atom at r0
c	gas cooling by molecular lines from Hollenbach & McKee (1979)
          colin = 1.0d+21
          colup = 1.0d+21
          xco = xco0
          xwat = xwat0
          cgas = cooling(dens,tgi,1,1,1)
          popsum = zero
          do ia = 1,na
            popa = popd(ia)*pi*ad(ia)**2
            popsum = popsum+popa
          enddo
          acgd = four*dens**2*popsum*sqrt(bk*tgi/xmg)*acom*2.5*bk
          hxr = pxr*axr/(fpi*r0**2)
          haccr = gmmdot/(fpi*r0**3*sigma0)
          print '(" gas heating and cooling at r0:",1p,4e10.2)',
     &      hxr,haccr,cgas,acgd
        endif

        if(doprint) print '(" calculate C&G disk model")'
        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)

c	calculate sublimation front and C&G disk structure
c	don''t know how to farm out
        ll = 1
        lw = 1
        ffstar = one
        nit = 0
c	rim location on the midplane
c	I think my approximations break down if qcol is too big
        if(qcol.lt.1.0d+01) then
          ginv = 0.4/gsubl-qcol  !  20-qcol
          powz = gsubl/0.2  !  0.1
          powr = two-qcol*powz  !  2-qcol/10
        else
          ginv = 0.4/gsubl-10.0  !  10
          powz = gsubl/0.2  !  0.1
          powr = two-10.0*powz  !  1
        endif
        ustar = four*tstar/hck
        iustar = slog(ustar/ud(1))/slog(ud(2)/ud(1))+1.5
        iustar = min(nu,max(1,iustar))
        ias1 = na/2+1
        ias = slog(asi/ad(ias1))/slog(ad(ias1+1)/ad(ias1))+1.5
        ias = min(na,max(ias1,ias))
        tmid = 1.0d+03
        hr2 = (tmid/tgrav)*rin/rstar
        hh2 = hr2*rin**2
        sigma = sigma0/((rin/r0)**pcol+(r0/rin)**qcol)
        pmid = tmid*sigma/sqrt(tpi*hh2)
        tsubli = tsubl*max(one,pmid)**gsubl
        usubl = four*tsubli/hck
        iusubl = slog(usubl/ud(1))/slog(ud(2)/ud(1))+1.5
        iusubl = min(nu,max(1,iusubl))
c	rim0: w/o disk radn,  rimi: current w/ disk radn, rim1: last w/
        rim0 = (rstar/two)*(teff/tsubli)**2
     &    *(phab(iustar,ias)/phab(iusubl,ias))**(one/powr)
        if(rim0.lt.rin) then
          if(doprint) print '("Warning: rim0 inside of rin")'
          rim0 = rin
        endif
        rim = rim0*(one+enhanc)**(one/powr)
        grmax = 0.99
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          if(igr.le.1) then
            ifrim = .false.
            zr1 = zero
            gr1 = one
            graz = one
            zrim = zero
            rphot = zero
            zphot = zero
            zrlast = zero
            do it = 1,nt
              do iu = 1,nu
                taud(iu,it) = zero
              enddo
              dtac4(it) = zero
            enddo
          endif
          do ir = kr,lr
            r = dr*(ir-half)
            if((r.gt.rin).and.(r.lt.rout)) then
              sigma = sigma0/((r/r0)**pcol+(r0/r)**qcol)
              tst4 = ostar*teff**4
              tir4 = graz*abss*tst4/two
              tir = tir4**quart
	tir = 1.0d+03
              hr2 = (tir/tgrav)*(r/rstar)
              hh2 = hr2*r**2
              pmid = tir*sigma/sqrt(tpi*hh2)
              tsubli = tsubl*max(one,pmid)**gsubl
            else
              sigma = zero
              pmid = zero
              tsubli = tsubl
            endif
c	find z rim semianalytically
            if(r.lt.rim) then
              if((pmid.gt.one).and.(ig.ne.1)) then
                usubl = four*tsubli/hck
                iusubl = slog(usubl/ud(1))/slog(ud(2)/ud(1))+1.5
                iusubl = min(nu,max(1,iusubl))
                rim0 = (rstar/two)*(teff/tsubli)**2
     &            *(phab(iustar,ias)/phab(iusubl,ias))**(one/powr)
                rim = rim0*(one+enhanc)**(one/powr)
              endif
              zri = zero
              gri = one
            else
              if(doprint.and.(.not.ifrim))
     &          print '(" found rim at",3f8.4,f8.2,1p,e10.2)',
     &            rim/au,rim0/au,rin/au,tsubli,pmid/tsubli
              ifrim = .true.
              alogr = slog(r/rim0)
              rr0 = (r/rim0)**powr
c	do I want to use gr1, or graz, here?
              arg = max(zero,ginv*alogr-slog(one+gr1*enhanc)/powz)
              if(zr1.eq.zero) then
c	first point beyond rim
                zri = sqrt(arg*hr2)
              else
c       z rim with last value of gr1
                zri = sqrt(arg*hr2)
c       solve for z rim iteratively
                do iit = 1,4
c       sin(graz)
                  gri = (rr0*exp(-powz*zri**2/hr2)-one)/enhanc
                  gri = min(grmax,gri)
c       tan(graz)
                  gri = gri/sqrt(one-gri**2)
                  zri = (zri+zr1+gri*dr/r)/two
                enddo
c       z rim including disk radiation
                gri = gri/sqrt(one+gri**2)
                arg = ginv*alogr-slog(one+gri*enhanc)/powz
                if(arg.gt.zero) zri = sqrt(arg*hr2)
              endif
              zrim = zri*r
            endif
            zr1 = zri
            gr1 = gri
            zphot = zero
            do it = 1,nt
              z = r*dcosth*(nt-it+half)
              if((rphot.gt.zero).and.(r.gt.rphot).and.(zphot.eq.zero)
     &          .and.(taud(iuv,it).gt.one)) then
                if((ig.eq.ng9).and.(it.eq.1)) then
                  zphot = r*min(coswind,zrlast*sqrt(slog(r/rphot)))
                elseif((ig.lt.ng9).and.(it.eq.itwind)) then
                  if(doprint) print '("Warning: dust opaque at thwind",
     &              3i6)', it,ir,ig
                  zphot = r*coswind
                else
                  dit = (taud(iuv,it)-one)/(taud(iuv,it)-taud(iuv,it-1))
                  zphot = r*dcosth*(nt-it+half+dit)
                  zrlast = (zphot/r)/sqrt(slog(r/rphot))
                endif
              endif
              if(z.gt.zrim) then
                condi = zero
              else
                condi = one
              endif
              do ia = 1,na
                condens(ia,it,ir,ig) = condi
              enddo
              opl(ml,it,ir,ig) = one
            enddo
            if((rphot.eq.zero).or.(r.lt.rphot).or.(zphot.eq.zero)) then
              graz = one
            else
              graz = zphot/(two*r*slog(r/rphot))
              graz = graz/sqrt(one+graz**2)
            endif
            if(r.le.rstar) then
              r = rstar
              ostar = one
            elseif(r.lt.(hun*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/two
            else
              ostar = quart*(rstar/r)**2
            endif
            tdrad = teff*ostar**fifth

c	C&G interior temperature
c	calculate conditions at midplane
            if(sigma.eq.zero) then
              tgi = tdrad
              tmid = tdrad
              pmid = zero
              dgmid = zero
              tst4 = (rstar/(two*r))**2*teff**4
              tir4 = graz*abss*tst4/two
            elseif(igr.eq.1) then
              tst4 = (rstar/(two*r))**2*teff**4
              tir4 = graz*abss*tst4/two
              tir = tir4**quart
	tir = 1.0d+03
              hr2 = (tir/tgrav)*(r/rstar)
              hh2 = hr2*r**2
              pmid = tir*sigma/sqrt(tpi*hh2)
              ir9 = ir*nn9-nn9/2
              tgi = (tg(nt,ir9,ng9)+tg(nt,ir9+1,ng9))/two
              tmid = tgi
              dgmid = pmid/tmid
              if(badnum(dgmid)) then
                if(doprint) print '("dgmid bad at",2i6,1p,6e9.2)',
     &            ir,ig,tmid,hh2,sigma,dgmid,tir,graz
                dgmid = zero
              endif
            else
c	radiation temperature of unattenuated stellar radiation
              tst4 = ostar*teff**4
c	divide by 2 b/c half of hot dust radiation goes back out
              tir4 = graz*abss*tst4/two
              tmid = (tir4+dtac4(nt))**quart
c	try leaving accr out of calculation of pmid
              tir = tir4**quart
	tir = 1.0d+03
              hr2 = (tir/tgrav)*(r/rstar)
              hh2 = hr2*r**2
              pmid = tir*sigma/sqrt(tpi*hh2)
              dgmid = pmid/tmid
              if(badnum(dgmid)) then
                print '("dgmid bad at",2i6,1p,6e10.2)',
     &            ir,ig,tir4,dtac4(nt),tmid,hh2,sigma,dgmid
                dgmid = zero
                pmid = zero
              endif
c	I think gastemp needs these
              dg(nt,ir,ig) = dgmid
              nit = -1
c	calculate gg including stellar radiation for gastemp
              do iu = 1,nu
                atten = sexp(-taud(iu,nt))
                tdhot = tsubli/two**quart
                gg(iu,nt,ir,ig) = bw(ud(iu),tmid)*(one-atten/two)
     &            +(atten/two)*bw(ud(iu),tdhot)*(tmid/tdhot)**4
     &            +ostar*atten*(bus(iu)+graz*rebus(iu))
              enddo
              taucont = taud(iuv,nt)
              if(ir.gt.kr) then
                colin = colr(nt,ir-1,ig)
                colup = colz(nt,ir-1,ig)
                tgi = tg(nt,ir-1,ig)
                tdi = td(nt,ir-1,ig)
              elseif(igr.gt.1) then
                colin = colr(nt,lr,ig-1)
                colup = colz(nt,lr,ig-1)
                tgi = tg(nt,lr,ig-1)
                tdi = td(nt,lr,ig-1)
              else
                colin = zero
                colup = zero
                tgi = tstar
                tdi = tstar
              endif
              if(taucont.gt.zero) then
                xwat = xwat0
                xco = xco0
              else
                xwat = xwat0/hun
                xco = xco0/hun
              endif
c	I probably should decrease vapor pressure if grains have settled
              pres = pmid
              tsubli = tsubl*max(one,pres)**gsubl
c	calculate midplane tg (and maybe sublimation)
              call gastemp(nt,ir,ig,iamax,tsubli,tgi,tdi)
              if(taucont.eq.zero) then
                tmid = min(tgi,tdrad)
              else
                tmid = min(tgi,tsubli)
              endif
              tmid = max(ten,tmid)
              if(r.lt.rin) pmid = zero
            endif  !  igr = 1?
c	iterate at each r
            nit = 0
  100       nit = nit-1
c	calculate hydro equil density, with temperature from one step down
            pres = pmid
c	tmid here is from last iteration
            dens = pmid/tmid
            ddfac = zero
            faccr = zero
c	step up from midplane
            do it = nt,1,-1
              costh = dcosth*(nt-it+half)
              z = r*costh
              dz = r*dcosth
              if((ig.eq.1).and.(it.gt.nt8)) then
c	average arrays from grid 9
                jr2 = ir*nn9
                jr1 = jr2-nn9+1
                jt2 = nt-(nt-it)*nn9
                jt1 = jt2-nn9+1
                sumdg = zero
                sumdd = zero
                sumtg = zero
                sumtd = zero
                sumph = zero
                sumopl = zero
                sumcol = zero
                do jt = jt1,jt2
                  do jr = jr1,jr2
                    sumdg = sumdg+dg(jt,jr,ng9)
                    sumdd = sumdd+dd(jt,jr,ng9)
                    sumtg = sumtg+tg(jt,jr,ng9)
                    sumtd = sumtd+td(jt,jr,ng9)
                  enddo
                  sumopl = sumopl+opl(ml,jt,jr2,ng9)
                  sumcol = sumcol+colr(jt,jr2,ng9)
                enddo
                dens = sumdg/xnn9
                ddfac = sumdd/xnn9
                tgi = sumtg/xnn9
                tdi = sumtd/xnn9
                pres = dens*tgi
                opl(ml,it,ir,ig) = sumopl/xn9
                colr(it,ir,ig) = sumcol/xn9
                if(sigma.gt.zero) then
                  faccr = faccr
     &              +gmmdot*dens*dz/(fpi*max(sigma,1.0d+18)*r**3)
                endif
                opapl = zero
                bwsum = zero
                do ia = 1,na
                  sumcon = zero
                  do jt = jt1,jt2
                    do jr = jr1,jr2
                      sumcon = sumcon+condens(ia,jt,jr,ng9)
                    enddo
                  enddo
                  condens(ia,it,ir,ig) = sumcon/xnn9
                enddo
                do iu = 1,nu
                  sumopa = zero
                  sumsca = zero
                  do jt = jt1,jt2
                    do jr = jr1,jr2
                      sumopa = sumopa+opa(iu,jt,jr,ng9)
                      sumsca = sumsca+sca(iu,jt,jr,ng9)
                    enddo
                  enddo
                  opai = sumopa/xnn9
                  opa(iu,it,ir,ig) = opai
                  scai = sumsca/xnn9
                  sca(iu,it,ir,ig) = scai
                  if(nit.eq.(-1))
     &              taud(iu,it) = taud(iu,it)+sqrt(opai*(opai+scai))*dr
                  bwdud = bw(ud(iu),tgi)*dud(iu)
                  if(opai.gt.zero) then
                    opapl = opapl+bwdud/opai
                  else
                    opapl = zero
                    goto 112
                  endif
                  bwsum = bwsum+bwdud
                enddo
                opapl = bwsum/opapl
                if(badnum(opapl)) then
                  if(doprint) print '("at",3i6," opapl =",1p,5e10.2)',
     &              it,ir,ig,opapl,bwsum,tgi,bwdud,opai
                  opapl = zero
                endif
  112           dtac4(it) = thrq*faccr*opapl*dz/sigsb
              elseif((sigma.eq.zero).or.(abs(costh).ge.coswind)) then
                dens = zero
                pres = zero
                tgi = tdrad
                tdi = tdrad
                ddfac = zero
                do iu = 1,nu
                  opa(iu,it,ir,ig) = zero
                  sca(iu,it,ir,ig) = zero
                  swd(iu,it,ir,ig) = ostar*bus(iu)
                  gg(iu,it,ir,ig) = ostar*bus(iu)
                  taud(iu,it) = zero
                enddo
                dtac4(it) = zero
                do ia = 1,na
                  condens(ia,it,ir,ig) = zero
                enddo
                opl(ml,it,ir,ig) = one
                colr(it,ir,ig) = zero
              else  !  dens.ne.0 and ig.ne.1
c	start with tgi from last iteration or neighboring cell
                if(nit.lt.(-1)) then
                  tgi = tg(it,ir,ig)
                  tdi = td(it,ir,ig)
                  taufac = opl(ml,it,ir,ig)
                elseif(it.lt.nt) then
                  tgi = tg(it+1,ir,ig)
                  tdi = td(it+1,ir,ig)
                  taufac = opl(ml,it+1,ir,ig)
                elseif(ir.gt.kr) then
                  tgi = tg(it,ir-1,ig)
                  tdi = td(it,ir-1,ig)
                  taufac = opl(ml,it,ir-1,ig)
                elseif(igr.gt.1) then
                  tgi = tg(it,lr,ig-1)
                  tdi = td(it,lr,ig-1)
                  taufac = opl(ml,it,lr,ig-1)
                else
                  tgi = tstar
                  tdi = tstar
                  taufac = one
                endif
c	maybe tgrav should change when H2 becomes H
                if((taufac.eq.one).or.(taud(iuv,it).gt.four)) then
                  hr2 = (tgi/tgrav)*(r/rstar)
                elseif(it.eq.nt) then
c	start with hr2 using tmid in thinned rim layer
                  hr2 = (tmid/tgrav)*(r/rstar)
                else
c	keep last good hr2?
                  continue
                endif
                hh2 = hr2*r**2
c	pressure change in half cell
                gaus = sexp(-(z*dz)/(two*hh2))
c	pressure at center of cell
                pres = pres*gaus
                if(pres.lt.one) pres = zero
c	density at center of cell
                dens = pres/tgi
c	dust density decrease exponential factor for a = 1 cm
                if(avisc*dens.gt.zero) then
                  ddfac = ddfac+z*dz/(avisc*hh2*sqrt(hh2)*dens)
                else
                  ddfac = zero
                endif
c	accumulated heat flux from viscous heating
                if(sigma.gt.zero) then
                  faccr = faccr
     &              +gmmdot*dens*dz/(fpi*max(sigma,1.0d+18)*r**3)
                endif
c	add accretion heating to heating by hot dust emission
c	tir4 is from midplane calculation
                tir = (tir4+dtac4(it))**quart
                tsubli = tsubl*max(pres,one)**gsubl
c	should I instead set condens = 0?
                tir = min(tir,tsubli)
                do iu = 1,nu
                  attenr = sexp(-taud(iu,it))
                  attenz = sexp(-taud(iu,it)*graz)
                  if(attenr*attenz.gt.zero) then
                    atten = one/(one/attenr+one/attenz)
                  else
                    atten = max(attenr,attenz)
                  endif
                  tdhot = tsubli/two**quart
                  gg(iu,it,ir,ig) = bw(ud(iu),tir)*(one-atten/two)
     &              +bw(ud(iu),tdhot)*(tir4/tdhot**4)*atten/two
     &              +ostar*(bus(iu)*attenr+graz*rebus(iu)*atten)
c	does this include downward scattered light?
                enddo
                taucont = taud(iuv,it)
                dg(it,ir,ig) = dens
                dd(it,ir,ig) = ddfac
                if(ir.gt.kr) then
                  colin = colr(it,ir-1,ig)
                  colup = colz(it,ir-1,ig)
                elseif(igr.gt.1) then
                  colin = colr(it,lr,ig-1)
                  colup = colz(it,lr,ig-1)
                else
                  colin = zero
                  colup = zero
                endif
c	should use better chemistry
                if(taucont.gt.zero) then
                  xwat = xwat0
                  xco = xco0
                else
                  xwat = xwat0/hun
                  xco = xco0/hun
                endif
c	calculate dust temperatures and opacities
                call gastemp(it,ir,ig,iamax,tsubli,tgi,tdi)
                dm(it,ir,ig) = iamax
                dens = pres/tgi
                if(dens.lt.tiny) then
                  dens = zero
                  pres = zero
                endif
c	calculate Planck mean opacity (should use Rosseland?)
                opapl = zero
                bwsum = zero
                do iu = 1,nu
                  opai = opa(iu,it,ir,ig)
                  if(opai.gt.zero) then
                    if(nit.eq.(-1)) taud(iu,it) = taud(iu,it)
     &                +sqrt(opai*(opai+sca(iu,it,ir,ig)))*dr
                    bwdud = bw(ud(iu),tgi)*dud(iu)
                    opapl = opapl+bwdud/opai
                    bwsum = bwsum+bwdud
                  else
                    opapl = zero
                    goto 110
                  endif
                enddo
                opapl = bwsum/opapl
                if(badnum(opapl)) then
                  if(doprint) print '("at",3i6," opapl =",1p,5e10.2)',
     &              it,ir,ig,opapl,bwsum,tgi,bwdud,opai
                  opapl = zero
                endif
  110           dtac4(it) = thrq*faccr*opapl*dz/sigsb
                taufac = opl(ml,it,ir,ig)
c	use recalculated tgi for hydro equil except at rim
                if((taufac.eq.one).or.(taud(iuv,it).gt.four)) then
                  hr2 = (tgi/tgrav)*(r/rstar)
                elseif(it.eq.nt) then
                  hr2 = (tmid/tgrav)*(r/rstar)
                endif
                hh2 = hr2*r**2
c	density change in half cell
                gaus = sexp(-(z*dz)/(two*hh2))
c	pressure at top of cell
                pres = pres*gaus
                if(pres.lt.one) pres = zero
              endif  !  ig.eq.1?
c	values at cell center
              dg(it,ir,ig) = dens
              dd(it,ir,ig) = ddfac
              tg(it,ir,ig) = tgi
              td(it,ir,ig) = tdi
              colr(it,ir,ig) = colin+taufac*dens*dr
              sumcon = zero
              do ia = 1,na
                sumcon = sumcon+condens(ia,it,ir,ig)*sexp(-ad(ia)*ddfac)
              enddo
c              dm(it,ir,ig) = sumcon*taufac
            enddo  !  it
            if((rphot.eq.zero).and.(taud(iuv,nt).gt.one)) then
              dtau = (opa(iuv,nt,ir,ig)+sca(iuv,nt,ir,ig))*dr
              dir = (taud(iuv,nt)-one)/dtau
              rphot = dr*(ir+half-dir)
            endif
c	save tmid for next iteration
            tmid = (9.*tg(nt,ir,ig)-tg(nt-1,ir,ig))/8.
c	integrate dtac4 and colz down through the disk
            dtaci = faccr/(two*sigsb)
            dtac4(1) = dtaci
            colup = zero
            colz(1,ir,ig) = zero
            do it = 2,nt
              dtaci = dtaci+dtac4(it)
              dtac4(it) = dtaci
              colup = colup+opl(ml,it,ir,ig)*dg(it,ir,ig)*dz
              colz(it,ir,ig) = colup
            enddo
            iw = ir+(ig-2)*nr/2
            if((ig.le.ng).and.(iw.ge.1).and.(iw.le.mw)) then
              smax = 1.0d+04
              spec(iw,4) = slog10(max(one,min(smax,tg(nt,ir,ig))))
              spec(iw,16) = spec(iw,4)
            endif
            if(nit.gt.-kit) goto 100
          enddo  !  ir
          if(ig.lt.ng) then
            tg(nt,nr/2,ig+1) = (tg(nt,nr-1,ig)+tg(nt,nr,ig))/two
          endif
        enddo  !  igr
        call copyarr(dm)
        call copyarr(dg)
        call copyarr(dd)
        call copywta(tg,dd)
        call copywta(td,dd)
        call copyarrr(opl,ml,ml)
        call copyarrr(opa,mu,nu)
        call copyarrr(sca,mu,nu)
        call copyarrr(swd,mu,nu)
        call copyarrr(gg,mu,nu)

        if(doprint) then
          print '(" midplane sublimation, photosphere radii:",
     &      2f10.4)', rim/au,rphot/au
          write(16,'(" midplane sublimation, photosphere radii:",
     &      2f10.4)') rim/au,rphot/au
        endif

c	plot mean intensity / stellar on xy grid
        if(doplot) then
c          print '(" plotting settling*condensation*taufactor")'
          print '(" plotting largest sublimating grain")'
          call plotcon(dm,.false.,'')
          call waitasec(wait,2)
          print '(" plotting dust heating radiation / stellar")'
          if(rim.lt.rgrid1) then
            call plotggfs(gg,ng9,iupl)
            call waitasec(wait,2)
          endif
          call plotggfs(gg,igpl,iupl)
          call waitasec(wait,2)
        endif
c	plot gas density and dust temperature
        if(doplot) then
          print '(" plotting log_10(hydrogen density)")'
          call plotcon(dg,.true.,'')
          call waitasec(wait,2)
          print '(" plotting C&G gas temperature")'
          call plotcon(tg,.false.,'')
        endif

        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)
        lw = nu
        sumfs = sumfs*astar+gmmdot/(two*rim*fpi*(dx0*two**ng)**2)

        nit = 0
        moltype = -abs(moltype)
        ffstar = one

c	loopback point for dust radiative transfer
  200   nit = nit+1

        if(doprint) print '(" calculate hydro equil disk structure")'
c	calculate hydro equil disk structure with new temperatures
        do it = 1,nt
          taud(iuv,it) = zero
        enddo
        rphot = zero
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
c	calculate disk structure
c	don''t know how to farm out
          do ir = kr,lr
            r = (ir-half)*dr
            taudr = taud(iuv,nt)
            if((taudr.gt.zero).and.(taudr.lt.four)) then
              irr = ir-1
  210         irr = irr+1
              taudr = taudr
     &          +(opa(iuv,it,irr,ig)+sca(iuv,it,irr,ig))*dr
c	take tmid from where taud = 4 along this los
c	so disk doesn''t expand in hot layer
              if((taudr.lt.four).and.(irr.lt.lr)) goto 210
              tmid = (9.*tg(nt,irr,ig)-tg(nt-1,irr,ig))/8.
            else
              tmid = (9.*tg(nt,ir,ig)-tg(nt-1,ir,ig))/8.
            endif
	tmid = 1.0d+03
            hh2 = (tmid/tgrav)*(r**3/rstar)
            hh0 = hh2
            if((r.gt.rin).and.(r.lt.rout)) then
              sigma = sigma0/((r/r0)**pcol+(r0/r)**qcol)
            else
              sigma = zero
            endif
            dgmid = sigma/sqrt(tpi*hh2)
            if(badnum(dgmid)) dgmid = zero
            dens = dgmid
            pres = dens*tmid
            ddfac = zero
c	step up from midplane
            do it = nt,1,-1
              costh = dcosth*(nt-it+half)
              if((ig.eq.1).and.(it.gt.nt8)) then
                jr2 = ir*nn9
                jr1 = jr2-nn9+1
                jt2 = nt-(nt-it)*nn9
                jt1 = jt2-nn9+1
                sumdg = zero
                sumdd = zero
                sumph = zero
                sumpd = zero
                do jr = jr1,jr2
                  do jt = jt1,jt2
                    dens = dg(jt,jr,ng9)
                    sumdg = sumdg+dens
                    sumdd = sumdd+dd(jt,jr,ng9)
                    sumph = sumph+dens*tg(jt,jr,ng9)
                  enddo
                enddo
                dens = sumdg/xnn9
                pres = sumph/xnn9
                ddfac = sumdd/xnn9
              elseif(abs(costh).ge.coswind) then
                dens = zero
                pres = zero
                ddfac = zero
              else
                z = r*costh
                dz = r*dcosth
                taudr = taud(iuv,it)
                if((opl(ml,it,ir,ig).lt.one).and.(taudr.lt.four)) then
                  irr = ir-1
  220             irr = irr+1
                  taudr = taudr
     &              +(opa(iuv,it,irr,ig)+sca(iuv,it,irr,ig))*dr
                  if((opl(ml,it,irr,ig).lt.one).and.(taudr.lt.four)
     &              .and.(irr.lt.lr)) goto 220
                  tgi = tg(it,irr,ig)
                  tdi = td(it,irr,ig)
                else
                  tgi = tg(it,ir,ig)
                  tdi = td(it,ir,ig)
                endif
c	tgrav should change when H2 becomes H
                hh2 = (tgi/tgrav)*(r**3/rstar)
c	dust density decrease exponential factor for a = 1 cm
                if(avisc*dens.gt.zero) then
                  ddfac = ddfac+z*dz/(avisc*hh2*sqrt(hh2)*dens)
                else
                  ddfac = zero
                endif
c	density change in half cell
                gaus = sexp(-(z*dz)/(two*hh2))
c	pressure at center of cell
                pres = pres*gaus
                if(pres.lt.one) pres = zero
                dens = pres/tgi
c	pressure at top of cell
                pres = pres*gaus
              endif  !  ig.eq.1?
              dg(it,ir,ig) = dens
              dm(it,ir,ig) = dens*(r/r0)**2
              dd(it,ir,ig) = ddfac
              taud(iuv,it) = taud(iuv,it)
     &          +(opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
            enddo  !  it
            colup = zero
            do it = 1,nt
              colup = colup+dg(it,ir,ig)*dz
              colz(it,ir,ig) = colup
            enddo
            if((rphot.eq.zero).and.(taud(iuv,nt).gt.one)) then
              dtau = (opa(iuv,nt,ir,ig)+sca(iuv,nt,ir,ig))*dr
              dir = (taud(iuv,nt)-one)/dtau
              rphot = dr*(ir+half-dir)
            endif
          enddo  !  ir
        enddo  !  igr
        call copyarr(dg)
        call copyarr(dm)
        call copyarr(dd)

c	plot gas density
        if(doplot) then
          print '(" plotting log_10(hydrogen density * r^2)")'
          monfile = ''
          if(nit.eq.(lit-1)) monfile = 'dens.con'
          call plotcon(dm,.true.,monfile)
          call waitasec(wait,2)
        endif

c	dust radiative transfer

        if(doprint) then
           print '(/" calculate dust radiative transfer  nit =",i3)',
     &       nit
           write(16,'(" calculate dust radiative transfer  nit =",i3)')
     &       nit
        endif

c	zero mean intensity arrays and spectrum
        call zerog(nu)
        do irot = 1,nrot
          do iu = 1,nu
            fw(iu,1,rot) = zero
          enddo
        enddo

c	loop over ray directions (grid orientations), farming out
      do irot0 = 1,nrot,mysize
        irot = irot0+myrank
        cz = one-dcrot*(irot-half)
        rot = acos(cz)
        sz = sin(rot)

        call zerof(nu)

c	call outer grid routine, which loops over rays and calls inner grids
c	stepd is called from grids to do dust radiative transfer
c	and add contributions to mean intensity array
        call grid0(fwj,ny/2)

c	add contribution to spectrum from +y and save in fw
        sumfwj = zero
        fwmn = two*fwj(1)*dud(1)/fwmx
        do iu = 1,nu
          fwiu = two*fwj(iu)
          sumfwj = sumfwj+fwiu*dud(iu)
          fw(iu,1,irot) = fwiu
          fwnorm = fwiu*dud(iu)/fwmx
          fwj(iu) = slog10(max(fwmn,fwnorm))
        enddo

c	plot images and spectrum
        if(doplot) then
          call movie(nc,1)
          call sleep(2)
          call movie(nc,2)
          call sleep(2)
          print '(" plotting log disk SED for irot =",i4)', irot
          jrot = (irot-1)/mysize+1
          call plotspec(fwj,mw,nu,1,jrot,0.,'')
        endif

c	broadcast calculated spectra and add contributions to mean intensity
        if(ifmpi) then
          do irank = 0,mysize-1
            call mpi_bcast(fw(1,1,irot0+irank),nu,
     &      MPI_DOUBLE_PRECISION,irank,MPI_COMM_WORLD,ierr)
          enddo
          do igr = 0,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            do ir = kr,lr
              call mpi_allreduce(MPI_IN_PLACE,fn(1,ir,ig),nt,
     &          MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
              call mpi_allreduce(MPI_IN_PLACE,ff(1,1,ir,ig),mu*nt,
     &          MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
            enddo
          enddo
        endif
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            do it = 1,nt
              gn(it,ir,ig) = gn(it,ir,ig)+fn(it,ir,ig)
              do iu = 1,nu
                gg(iu,it,ir,ig) = gg(iu,it,ir,ig)+ff(iu,it,ir,ig)
              enddo
            enddo
          enddo
        enddo

c	end loop over grid orientations
      enddo  !  irot

c	fill in inner halves so I can smooth radially
        call copyarr(gn)
        call copyarrr(gg,mu,nu)
c	divide gg and qq by gn and average if appropriate
        call ggavg(nu,ntsm,nrsm)
        call copyarrr(gg,mu,nu)
c	plot gg before adding starlight
        if(doplot) then
          print '(" plotting dust heating emission / stellar")'
          call plotggfs(gg,ng9,iupl)
          call waitasec(wait,2)
        endif

c	radiative transfer of starlight (w/o sources)
c	note: not doing star for nt < irot < nnt
        if(doprint) print '(" starlight radiative transfer")'
        call zerof(nu)
        dxstar = rstar/four
        idx = 4
        ix1 = nx9/2+1-idx
        ix2 = nx9/2+idx
c	set iy to a value unused by dust rad xfer
        iy = nx9/2+1
        nrots = nt
c	why did I do an additional mysize rots?
        if(nt+mysize.le.nnt) nrots = nt+mysize
        do irot0 = 1,nrots,mysize
          irot = irot0+myrank
          cz = dcost*(nt-irot+half)
          rot = acos(cz)
          sz = sin(rot)
          do iu = 1,nu
            fw(iu,2,irot) = zero
          enddo
          sumbg = zero
          do ix = ix1,ix2
            x = dxstar*(ix-nx9/2-half)
            gwt = sqrt(one-(x/rstar)**2)
            sumbg = sumbg+bstar*gwt
c	start with stellar double blackbody spectrum
            do iu = 1,nu
              fwj(iu) = bus(iu)
            enddo
            do igr = 1,ng
              call startgr(igr,ig,kr,lr,lt,dr,dcosth)
              if(igr.eq.1) then
                lz = nz9
                dz = dz9
                kz = lz/2+1+nint(rstar/dz)
              else
                lz = nz
                kz = 3*lz/4+1
                twig = two**ig
                dz = dz0*twig
              endif
              do iz = kz,lz
                call stepd(fwj,ix,iy,iz,ig)
              enddo
            enddo  !  igr
c	sum stellar spectra in fw(iu,2,irot)
            do iu = 1,nu
              fw(iu,2,irot) = fw(iu,2,irot)+bstar*gwt*fwj(iu)
            enddo
          enddo  !  ix
          if(ifmpi) then
            do irank = 0,mysize-1
              call mpi_bcast(fw(1,2,irot0+irank),nu,
     &          MPI_DOUBLE_PRECISION,irank,MPI_COMM_WORLD,ierr)
            enddo
          endif
        enddo  !  irot

c	do finer rot steps for grid 9
        ig = ng9
        do ir = 1,nr9
          do it = 1,nt
            do iu = 1,nu
              ff(iu,it,ir,ig) = zero
            enddo
            fn(it,ir,ig) = zero
          enddo
        enddo
        igr = 0
        call startgr(igr,ig,kr,lr,lt,dr,dcosth)
        lz = nz9
        dz = dz9
        kz = lz/2+nint(rstar/dz+one)
        do irot0 = 1,nt,mysize
          irot = irot0+myrank
          cz = dcos9*(nt-irot+half)
          rot = acos(cz)
          sz = sin(rot)
          do ix = ix1,ix2
c	start with stellar double blackbody spectrum
            do iu = 1,nu
              fwj(iu) = bus(iu)
            enddo
            do iz = kz,lz
              call stepd(fwj,ix,iy,iz,ig)
            enddo
          enddo  !  ix
        enddo  !  irot

c	sum starlight mean intensities
        if(ifmpi) then
          do igr = 0,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            do ir = kr,lr
              call mpi_allreduce(MPI_IN_PLACE,fn(1,ir,ig),nt,
     &          MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
              call mpi_allreduce(MPI_IN_PLACE,ff(1,1,ir,ig),mu*nt,
     &          MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
            enddo
          enddo
        endif
c	add starlight to mean intensities
c	note: averaging over half of star at it = 1
        nzero = 0
        sumfn = zero
        sumn = zero
        fnmin = one
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            if(r.lt.rstar) then
              ostar = one
            elseif(r.lt.(hun*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/two
            else
              ostar = quart*(rstar/r)**2
            endif
            do it = 1,nt
              if(smooth.and.(igr.eq.0)) then
c	3 point angular average in grid 9
c	is this a good idea?
                jt1 = it-1
                jt2 = it+1
                fni = zero
                do jt = jt1,jt2
                  if(jt.lt.1) then
                    kt = 1-jt
                  elseif(jt.gt.nt) then
                    kt = 2*nt-jt+1
                  else
                    kt = jt
                  endif
                  fni = fni+fn(kt,ir,ig)
                enddo
              else
                fni = fn(it,ir,ig)
              endif
              if(fni.gt.fnmin) then
                sumfn = sumfn+fni
                sumn = sumn+one
                wtf = ostar/fni
                if(smooth.and.(igr.eq.0)) then
c	3 point angular smooth of ff in grid 9
                  jt1 = it-1
                  jt2 = it+1
                  do jt = jt1,jt2
                    if(jt.lt.1) then
                      kt = 1
                    elseif(jt.gt.nt) then
                      kt = 2*nt-jt+1
                    else
                      kt = jt
                    endif
c	put ffstar correction into gg and qq while adding starlight
                    do iu = 1,nu
                      gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/ffstar
     &                  +wtf*ff(iu,kt,ir,ig)
                    enddo
                  enddo
                else
                  do iu = 1,nu
                    gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/ffstar
     &                  +wtf*ff(iu,it,ir,ig)
                    ff(iu,it,ir,ig) = wtf*ff(iu,it,ir,ig)
                  enddo
                endif
              elseif((it.le.lt).and.(r.gt.rin)) then
                nzero = nzero+1
              endif
            enddo  !  it
          enddo  !  ir
        enddo  !  igr
        if(doprint)
     &    print '(" for starlight, mean rays/cell, zeros:",f8.2,i8)',
     &      sumfn/sumn,nzero

c	plot mean intensity / stellar on xy grid
        call copyarrr(gg,mu,nu)
        if(doplot) then
          if(rim.lt.rgrid1) then
            print '(" plotting dust heating starlight / unattentuated")'
            call plotggfs(ff,ng9,iupl)
            call waitasec(wait,2)
            call plotggfs(gg,ng9,iupl)
            call waitasec(wait,2)
          endif
          print '(" plotting dust heating radiation / stellar")'
          call plotggfs(gg,igpl,iupl)
          call waitasec(wait,2)
        endif

        if(doprint)
     &    print '(" calculate gas and dust temperatures")'
c	first follow radial rays in grid 9
c	loop over ray directions, farming out
        dtsum = zero
        dtmax = zero
        do it0 = 1,nt,mysize
          it = it0+myrank
          colin = zero
          taucont = zero
          igr = 0
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            dens = dg(it,ir,ig)
            if((r.gt.rin).and.(r.lt.rout)) then
              sigma = sigma0/((r/r0)**pcol+(r0/r)**qcol)
            else
              sigma = zero
            endif
            tgi = tg(it,ir,ig)
            tdi = td(it,ir,ig)
            tsubli = tsubl*max(one,pres)**gsubl
            colup = colz(it,ir,ig)
            if(taucont.eq.zero) then
              xwat = xwat0/hun
              xco = xco0/hun
            else
              xwat = xwat0
              xco = xco0
            endif
c	gastemp calculates temperatures of grain species and gas
c	from equations of heating and cooling
            call gastemp(it,ir,ig,iamax,tsubli,tgi,tdi)
            dti = tg(it,ir,ig)-tgi
            tg(it,ir,ig) = tgi
            td(it,ir,ig) = tdi
            colin = colin+dens*dr
            colr(it,ir,ig) = colin
            taucont = taucont+(opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
            tauv(it,ir,ig) = taucont
            dtsum = dtsum+dti**2
            if(abs(dti).gt.abs(dtmax)) dtmax = dti
          enddo  !  ir
c	broadcast after doing ir loop
          if(ifmpi) then
            do ir = kr,lr
              call mpi_allgather(
     &          tg(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          tg(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          td(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          td(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          colr(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          colr(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          tauv(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          tauv(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          opa(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &          opa(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          sca(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &          sca(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          swd(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &          swd(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
              call mpi_allgather(
     &          opl(1,it,ir,ig),ml,MPI_DOUBLE_PRECISION,
     &          opl(1,it0,ir,ig),ml,MPI_DOUBLE_PRECISION,
     &          MPI_COMM_WORLD,ierr)
            enddo  !  ir
          endif
        enddo  !  it0

c	follow or continue to follow rays on grids 1-8
c	loop over ray directions, farming out
        do it0 = 1,nt,mysize
          it = it0+myrank
          lt9 = nt-nn9*(nt-it)
          kt9 = lt9-nn9+1
          colin = zero
          taucont = zero
          do igr = 1,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            do ir = kr,lr
              r = dr*(ir-half)
              if(r.lt.rin) taucont = zero
              dens = dg(it,ir,ig)
              if((r.gt.rin).and.(r.lt.rout)) then
                sigma = sigma0/((r/r0)**pcol+(r0/r)**qcol)
              else
                sigma = zero
              endif
              tgi = tg(it,ir,ig)
              tdi = td(it,ir,ig)
              tsubli = tsubl*max(one,pres)**gsubl
              colup = colz(it,ir,ig)
              if(taucont.eq.zero) then
                xwat = xwat0/hun
                xco = xco0/hun
              else
                xwat = xwat0
                xco = xco0
              endif
              call gastemp(it,ir,ig,iamax,tsubli,tgi,tdi)
              if((ig.eq.ng9).or.((ig.gt.1).and.(it.gt.itwind))) then
                dti = tg(it,ir,ig)-tgi
                dtsum = dtsum+dti**2
                if(abs(dti).gt.abs(dtmax)) dtmax = dti
              endif
              tg(it,ir,ig) = tgi
              td(it,ir,ig) = tdi
              colin = colin+dens*dr
              colr(it,ir,ig) = colin
              taucont = taucont
     &          +(opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
              tauv(it,ir,ig) = taucont
            enddo  !  ir
            if((igr.eq.1).and.(kt9.ge.1).and.(lt9.le.nt)) then
              do ir = kr,lr
                lr9 = ir*nn9
                kr9 = lr9-nn9+1
                sumtg = zero
                sumtd = zero
                sumcolr = zero
                sumtauv = zero
                do it9 = kt9,lt9
                  do ir9 = kr9,lr9
                    sumtg = sumtg+tg(it9,ir9,ng9)
                    sumtd = sumtd+td(it9,ir9,ng9)
                  enddo
                  sumcolr = sumcolr+colr(it9,lr9,ng9)
                  sumtauv = sumtauv+tauv(it9,lr9,ng9)
                enddo
                tg(it,ir,ig) = sumtg/xnn9
                td(it,ir,ig) = sumtd/xnn9
                colin = sumcolr/xn9
                colr(it,ir,ig) = colin
                taucont = sumtauv/xn9
                tauv(it,ir,ig) = taucont
              enddo
            endif
          enddo  !  igr
          if(ifmpi) then
            do igr = 0,ng
              call startgr(igr,ig,kr,lr,lt,dr,dcosth)
              do ir = kr,lr
                call mpi_allgather(
     &            tg(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            tg(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            td(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            td(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            colr(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            colr(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            tauv(it,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            tauv(it0,ir,ig),1,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            opa(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &            opa(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            sca(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &            sca(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            swd(1,it,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &            swd(1,it0,ir,ig),mu,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
                call mpi_allgather(
     &            opl(1,it,ir,ig),ml,MPI_DOUBLE_PRECISION,
     &            opl(1,it0,ir,ig),ml,MPI_DOUBLE_PRECISION,
     &            MPI_COMM_WORLD,ierr)
              enddo  !  ir
            enddo  !  igr
          endif  !  mpi
        enddo  !  it

        call copyarr(tg)
        call copyarr(td)
        call copytau(colr)
        call copytau(tauv)
        call copyarrr(opa,mu,nu)
        call copyarrr(sca,mu,nu)
        call copyarrr(swd,mu,nu)
        call copyarrr(ff,mu,nu)
        call copyarrr(gg,mu,nu)
        call copyarrr(opl,ml,ml)

c	add starlight to dust spectra and compare outgoing to stellar
        ffstar = zero
        fnstar = zero
        starsum = zero
        nnt = nint(one/dcost)
        xtrot = real(nnt)/real(nrot)
        do irot = 1,nrot
          if(nnt.eq.nrot) then
            it = irot-nrot+nt
          else
            it = nint((irot-nrot-half)*xtrot+half)+nt
          endif
          sumfw = zero
          starf = zero
          if((it.ge.1).and.(it.le.nt)) then
            do iu = 1,nu
              fw(iu,1,irot) = fw(iu,1,irot)+fw(iu,2,it)
              sumfw = sumfw+fw(iu,1,irot)*dud(iu)
              starf = starf+fw(iu,2,it)*dud(iu)
            enddo
          else
            do iu = 1,nu
              fw(iu,1,irot) = fw(iu,1,irot)+astar*bus(iu)
              sumfw = sumfw+fw(iu,1,irot)*dud(iu)
              starf = starf+astar*bus(iu)*dud(iu)
            enddo
          endif
          ffstar = ffstar+sumfw
          starsum = starsum+starf
          fnstar = fnstar+one
        enddo
        ffstar = ffstar/(fnstar*sumfs)
        starsum = starsum/(fnstar*sumfs)
        if(doprint) then
          print '(" outgoing flux/stellar:",f10.3)', ffstar
          print '(" just starlight:",f10.3)', starsum
          write(16,'(" outgoing flux/stellar:",f10.3)') ffstar
          write(16,'(" just starlight:",f10.3)') starsum
        endif
        ffstar = ((ffstar+one)/two-starsum)/(one-starsum)
        if(nit.lt.8) then
          ffstar = one
        else
          ffstar = min(1.7,max(0.6,ffstar))
        endif
        if(doprint) print '(" correction factor:",f10.3)', one/ffstar

c	make plot of density and temperature
        nwspec = min(mw,ng*nr/2)
        do igr = 1,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = kr,lr
            r = dr*(ir-half)
            it = nt
            iw = ir+(ig-2)*nr/2
            if((ig.le.ng).and.(iw.gt.0).and.(iw.le.mw)) then
              if(r.lt.rstar) r = rstar
              dens = dg(it,ir,ig)
              tgi = tg(it,ir,ig)
              tdi = td(it,ir,ig)
              tsubli = tsubl*max(one,dens*tgi)**gsubl
              smin = one
              smax = 1.0d+08
              spec(iw,1) = slog10(max(smin,min(smax,
     &          dens**quart)))
              spec(iw,2) = slog10(max(smin,min(smax,tgi)))
              spec(iw,3) = slog10(max(smin,min(smax,tdi)))
              spec(iw,4) = slog10(max(smin,min(smax,tsubli)))
              spec(iw,5) = spec(iw,16)
            endif
          enddo  !  ir
        enddo  !  igr

        if(ifmpi) call mpi_allreduce(MPI_IN_PLACE,dtsum,1,
     &        MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
        rmsdt = sqrt(dtsum/(nr9*nt+(ng-1)*nr2*(nt-itwind)))
        if(doprint) then
          print '(" midplane sublimation, photosphere radii:",
     &      2f10.4)', rim/au,rphot/au
          write(16,'(" midplane sublimation, photosphere radii:",
     &      2f10.4)') rim/au,rphot/au
          print '(" rms, max dt:",1p,2e10.2)', rmsdt,dtmax
          write(16,'(" rms, max dt:",1p,2e10.3)') rmsdt,dtmax
        endif

        if(doplot) then
          print '(" plotting log(n_H, T_g, T_subl, T_C&G) vs. ir")'
          call plotspec(spec,mw,nwspec,5,1,0.,'')
          call sleep(4)
c	plot tauv
          print '(" plotting log tauv")'
          monfile = ''
          call plotcon(tauv,.true.,monfile)
          call waitasec(wait,2)
c	plot temperature
          print '(" plotting gas temperature")'
          monfile = ''
          if(nit.eq.lit) monfile = 'temp.con'
          call plotcon(tg,.false.,monfile)
          call waitasec(wait,2)
c	plot SEDs including starlight
          print '(" plotting SEDs including starlight")'
          fwmn = fw(1,1,1)*dud(1)/fwmx
          do ipl = 1,8
            irot = (ipl*nrot)/8
            do iu = 1,nu
              fwnorm = fw(iu,1,irot)*dud(iu)/fwmx
              spec(iu,ipl) = slog10(max(fwmn,fwnorm))
            enddo
            xn13 =-slog10((ud(36)*fw(36,1,irot))/(ud(30)*fw(30,1,irot)))
     &            /slog10(ud(36)/ud(30))
            xn25 =-slog10((ud(30)*fw(30,1,irot))/(ud(26)*fw(26,1,irot)))
     &            /slog10(ud(30)/ud(26))
            xn31 =-slog10((ud(30)*fw(30,1,irot))/(ud(24)*fw(24,1,irot)))
     &            /slog10(ud(30)/ud(24))
            print '(" SED indices:",3f8.3)', xn13,xn25,xn31
            write(16,'(" SED indices:",3f8.3)') xn13,xn25,xn31
          enddo
          if(nit.eq.lit) monfile = 'seds.spec'
          call plotspec(spec,mw,nu,8,1,0.,monfile)
          call waitasec(wait,2)
        endif
        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)

c	loop back to hydro equil calculation
        if(nit.lt.lit) goto 200

c	save results of dust calculation
        if(doprint) then
          call savemodel
          open(unit=10,file='seds.mon',access='sequential')
          do iu = 1,nu
            write(10,'(i6,8f9.5)') iu,(spec(iu,ipl),ipl=1,8)
          enddo
          close(unit=10)
          print 
     &      '(" midplane dissociation radius for H2 and CO:",2f10.3)',
     &      rhi/au,rphoto/au
          write(16,
     &      '(" midplane dissociation radius for H2 and CO:",2f10.3)')
     &      rhi/au,rphoto/au
        endif
        if(ifmpi) call mpi_barrier(MPI_COMM_WORLD,ierr)

        if(ifmpi) call mpi_finalize
        stop

      end



      subroutine startgr(igr,ig,kr,lr,lt,dr,dcosth)

        implicit real*8 (a-h,o-z)
        parameter (two = 2.0d+00)
        save dr0,dr9,dcost,dcos9
        save nr,nr3,nr9,ng9,nt,nt8

        if((igr.eq.0).or.(igr.eq.ng9)) then
          ig = ng9
          dr = dr9
          kr = 1
          lr = nr9
          dcosth = dcos9
        else
          ig = igr
          dr = dr0*two**ig
          kr = nr3
          lr = nr
          dcosth = dcost
        endif
        if(ig.eq.1) then
          kr = 1
          lt = nt8
        else
          lt = nt
        endif
        return

      entry initgr(xr0,xr9,xcost,xcos9,ir,ir9,ig9,it,it8)

        dr0 = xr0
        dr9 = xr9
        dcost = xcost
        dcos9 = xcos9
        nr = ir
        nr9 = ir9
        ng9 = ig9
        nr3 = nr/2+1
        nt = it
        nt8 = it8
        return

      end



      subroutine dustinit(xc,xsi,ac,asi,tuv,xuv,emiss,tsubli)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)
        common /taus/ taucont,tauline,tauterm,iuv,iusil,iufl,iwterm
        dimension spec(mu,ma)
        logical doprint,badnum
        parameter (pi = 3.14159265d+00)
        parameter (tpi = 6.28318530d+00)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (ten = 1.0d+01)
        parameter (hun = 1.0d+02)

        nr2 = nr/2
        doprint = (myrank.eq.mysize/2)
c	a = 0.3 nm - 100 um (may be changed)
c	na = 24; 12 each for carbon and silicates
        la = na/2
        arat = ten**half
        a0 = 1.0d-07
c	roll off distributions at ac and asi
c	grain distribution power laws, dn/dlna ~ a^-pow
        powc = 2.5d+00
        pows = 2.0d+00
c       assume AC grain density = 2 g/cc or 1e23 C atoms/cc
c       so grain vol/C = 1e-23, or vol/H = 1e-23*XC
c       solar XC = 3.6e-04, gas particles / H atom = 0.6
        xc = 3.6d-04*xc
        cvolg = (1.0d-23)*xc/(0.6d+00)
c	assume Sil grain density = 3.6 g/cc, with 172 AMU/Si,
c	or 1.25e22 Si/cc, or 8e-23 cc/Si
        xsi = 3.5d-05*xsi
        svolg = (8.0d-23)*xsi/(0.6d+00)
c	u = 11.5 - 10^5 cm-1
        u0 = ten
        dlogu = 4.0d+00/nu
        urat = ten**dlogu
        durat = (urat-one)/sqrt(urat)
c	w = 328 - 3162 cm-1
        w0 = ten**2.5
        dlogw = one/nu
        wrat = ten**dlogw
        dwrat = (wrat-one)/sqrt(wrat)
        wv = 2.0d+04
        iuv = slog10(wv/u0)/dlogu+half
        wfl = 5.0d+04
        iufl = slog10(wfl/u0)/dlogu+half
        wsil = 1.0d+03
        iusil = slog10(wsil/u0)/dlogu+half
        iwsil = slog10(wsil/w0)/dlogw+half
        if(iwsil.ne.iusil)
     &    print '("Warning: iwsil.ne.iusil",2i6)',iwsil,iusil
        sumv = zero
        if(na.ne.24) print '("Warning: na =",i4)', na
        volc = zero
        vols = zero
        a = a0
        do ia = 1,la
          iac = ia
          ias = ia+la
          a = a*arat
          ad(iac) = a
          ad(ias) = a
c	approximation to Weingartner & Draine (2001)
c	note: popd = dn/dlna = a*dn/da since a is spaced logarithmically
          pops = one/((a/asi)**pows+(a/asi)**10)
          popc = one/((a/ac)**powc+(a/ac)**10)
          popd(iac) = popc
          popd(ias) = pops
          vola = (4.*pi/3.)*a**3
          volc = volc+popc*vola
          vols = vols+pops*vola
        enddo
        cnorm = cvolg/volc
        snorm = svolg/vols
        do ia = 1,la
          iac = ia
          ias = ia+la
          popd(iac) = cnorm*popd(iac)
          popd(ias) = snorm*popd(ias)
        enddo
        u = u0
        do iu = 1,nu
          u = u*urat
          ud(iu) = u
          dud(iu) = u*durat
          bus(iu) = bw(u,tstar)+xuv*(tstar/tuv)**4*bw(u,tuv)
        enddo
c       initialize phab, phsc, reading from ac.out and as.out
        open(unit=10,file='ac.out',access='sequential')
        read(10,'(1x)')
        do ia = 1,la
          a = ad(ia)
c	skip every other a
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
          do iu = nu,1,-1
c	read every 4th line
            read(10,'(6e12.4///)') amie,wlen,qext,qabs,qscat,cosscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one)
     &        .or.(abs(1.0e+04-ud(iu)*wlen).gt.one)))
     &        print '(2i6,1p,4e10.2)', ia,iu,a,amie,ud(iu),wlen
            phab(iu,ia) = qabs
            if((cosscat.gt.zero).and.(cosscat.lt.one))
     &        qscat = (one-cosscat)*qscat
            phsc(iu,ia) = qscat
          enddo
        enddo
        close(unit=10)
        open(unit=10,file='as.out',access='sequential')
        read(10,'(1x)')
        do ia = la+1,na
          a = ad(ia)
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
          do iu = nu,1,-1
            read(10,'(5e12.4///)') amie,wlen,qext,qabs,qscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one)
     &        .or.(abs(1.0e+04-ud(iu)*wlen).gt.one)))
     &        print '(2i6,1p,4e10.2)', ia,iu,a,amie,ud(iu),wlen
            phab(iu,ia) = qabs
            phsc(iu,ia) = qscat
          enddo
        enddo
        close(unit=10)
c       finer grid for last iteration
        w = w0
        do iu = 1,nu
          w = w*wrat
          wd(iu) = w
        enddo
        open(unit=10,file='ac.out',access='sequential')
        read(10,'(1x)')
        do ia = 1,la
          a = ad(ia)
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
c	skip to 31.6um
          do iu = 1,96
            read(10,'(1x)')
          enddo
          do iu = nu,1,-1
            read(10,'(5e12.4)') amie,wlen,qext,qabs,qscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one)
     &        .or.(abs(1.0e+04-wd(iu)*wlen).gt.one)))
     &        print '(2i6,1p,4e10.2)', ia,iu,a,amie,wd(iu),wlen
            phid(iu,ia) = qabs
          enddo
          do iu = 1,96
            read(10,'(1x)')
          enddo
        enddo
        close(unit=10)
        open(unit=10,file='as.out',access='sequential')
        read(10,'(1x)')
        do ia = la+1,na
          a = ad(ia)
          do iu = nu,1,-1
            read(10,'(///)')
          enddo
          do iu = 1,96
            read(10,'(1x)')
          enddo
          do iu = nu,1,-1
            read(10,'(5e12.4)') amie,wlen,qext,qabs,qscat
            if(doprint.and.((abs(1.0e+04-amie/a).gt.one)
     &        .or.(abs(1.0e+04-wd(iu)*wlen).gt.one)))
     &        print '(2i6,1p,4e10.2)', ia,iu,a,amie,wd(iu),wlen
            phid(iu,ia) = qabs
          enddo
          do iu = 1,96
            read(10,'(1x)')
          enddo
        enddo
        close(unit=10)

c	calculate emis = 1-albedo averaged over stellar spectrum
c	but note that sqrt(avg) != avg(sqrt)
        emissum = zero
        emissil = zero
        sumbus = zero
        do iu = 1,nu
          busdud = bus(iu)*dud(iu)
          opac = zero
          scat = zero
          opacsil = zero
          scatsil = zero
          do ia = 1,na
            a = ad(ia)
            popsig = popd(ia)*pi*a**2*phab(iu,ia)
            if(mod(ia,2).eq.0) spec(iu,ia/2) = slog10(popsig)+2.4d+01
            opac = opac+popsig
            if(ia.gt.la) opacsil = opacsil+popsig
            popsig = popd(ia)*pi*a**2*phsc(iu,ia)
            scat = scat+popsig
            if(ia.gt.la) scatsil = scatsil+popsig
          enddo
          emis(iu) = opac/(opac+scat)
          emissum = emissum+busdud*opac/(opac+scat)
          emissil = emissil+busdud*opacsil/(opacsil+scatsil)
          sumbus = sumbus+busdud
        enddo
        emiss = emissum/sumbus
        emissil = emissil/sumbus
        if(doprint) print '(" albedos and reflectivities:",4f8.3)',
     &    one-emiss,one-emissil,one-sqrt(emiss),one-sqrt(emissil)
        treem = tsubli/two**quart
        reemiss = half*sqrt(emiss)*(treem/tstar)**4
        do iu = 1,nu
          rebus(iu) = bus(iu)
     &      +sqrt(emis(iu))*(reemiss*bw(ud(iu),treem)-bus(iu))
        enddo

c	plot opacity and populations
        if(doprint) then
          call plotspec(spec,mu,nu,na/2,1,0.,'')
          call sleep(4)
          do iu = 1,nu
            do ia = 2,na,2
              spec(iu,ia/2) = max(-5.,slog10(phab(iu,ia)))
            enddo
          enddo
          call plotspec(spec,mu,nu,na/2,1,0.,'')
        endif

        return
      end



      subroutine stepd(fwi,ix,iy,iz,igr)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)

        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /xj/ xju,xjp,xjq,xjr,xk
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)
        common /taus/ taucont,tauline,tauterm,iuv,iusil,iufl,iwterm
        common /exptab/ expi(10239),too
        dimension fwi(mw),ffi(ml)
        logical good,doprint,badnum
        logical starlt

        parameter (pi = 3.14159265d+00)
        parameter (sqrtpi = 1.772454d+00)
        parameter (sqrt2 = 1.414214d+00)
        parameter (alog2 = 0.6931472d+00)
        parameter (small = 1.0d-03)
        parameter (tiny = 1.0d-06)
        parameter (zilch = 1.0d-200)
        parameter (hun = 1.0d+02)
        parameter (big = 5.0d+03)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (eighth = 0.125d+00)
        parameter (au = 1.496d+13)

        starlt = (iy.gt.nx9/2)
        jg = abs(igr)
        if(starlt) then
          lx = nx9
          dx = rstar/four
          if((jg.eq.1).or.(jg.eq.ng9)) then
            jg = ng9
            lz = nz9
            dz = dz9
          else
            twig = two**jg
            lz = nz
            dz = dz0*twig
          endif
        elseif((jg.eq.1).or.(jg.eq.ng9)) then
          jg = ng9
          lx = nx9
          dx = dx9
          lz = nz9
          dz = dz9
        else
          twig = two**jg
          lx = nx
          dx = dx0*twig
          lz = nz
          dz = dz0*twig
        endif
        ly = lx/2
c	finer steps in middle of grid
        if(igr.lt.0) dx = dx/two
        doprint = (myrank.eq.0).and.(ix.eq.48).and.(iy.eq.ly)
     &    .and.(jg.eq.ng9)
        dx2 = dx/two
        dz2 = dz/two
        ds = dz
        x = dx2*(2*ix-lx-1)
        y = dx2*(2*iy-lx-1)
        z = dz2*(2*iz-lz-1)
        r = sqrt(x**2+y**2+z**2)
        cz = cos(rot)
        sz = sin(rot)
        costh = abs(cz*z+sz*x)/r
        if(costh.gt.one)
     &    print '(" In stepd x,z,sz,cz,costh =",1p,5e10.2)',
     &      x,z,sz,cz,costh
        sinth = sqrt(one-costh**2)
        dcosth = dcost
        cosmx = nt*dcos9
        ig = slog(r/(nr*dr0))/alog2+one
        if(jg.le.ng) then
          ig = max(ig,jg)
          if(ig.eq.jg) then
            dr = dr0*twig
            lr = nr
            gwt = one
          elseif((ig.eq.jg+1).and.(ig.le.ng)) then
            dr = two*dr0*twig
            lr = nr
            gwt = eighth
          else
            return
          endif
        else
          ig = max(ig,1)
          if(ig.gt.1) then
            ig = 2
            dr = four*dr0
            lr = nr
            gwt = (dz9/(four*dz0))*(dx9/(four*dx0))**2
          elseif(abs(costh).gt.cosmx) then
            ig = 1
            dr = two*dr0
            lr = nr
            gwt = (dz9/(two*dz0))*(dx9/(two*dx0))**2
          else
            ig = ng9
            dr = dr9
            lr = nr9
            gwt = one/nn9
            dcosth = dcos9
          endif
        endif
        ri = r/dr+half
        ir = nint(ri)
        if((ir.lt.1).or.(ir.gt.lr)) return
        dir = ri-ir
        if(dir.lt.zero) then
          dir = -dir
          jr = max(ir-1,1)
        else
          jr = min(ir+1,lr)
        endif
        ti = nt-costh/dcosth+half
        it = nint(ti)
        if((it.lt.1).or.(it.gt.nt)) return
        dit = ti-it
        if(dit.lt.zero) then
          dit = -dit
          jt = max(it-1,1)
        else
          jt = min(it+1,nt)
        endif
        dirt = one-dir-dit
        if(igr.lt.0) gwt = gwt/four
        if(iy.gt.ly) then
          if(abs(x).lt.rstar) then
            gwt = sqrt(one-(x/rstar)**2)
          else
            gwt = zero
          endif
          if(r.lt.rstar) return
        endif

c	let star absorb, but not emit
c	put starlight in at start of starlight rays
        if(starlt) then
          continue
        elseif((jg.eq.ng9).and.(z.gt.zero).and.(z.lt.dz)) then
c	crossing midplane
          taucont = zero
          xi = abs(x)-dx2
          yi = abs(y)-dx2
          ri = sqrt(xi**2+yi**2)
          if(ri.lt.rstar)then
c	inner corner passes through star
            xo = xi+dx
            yo = yi+dx
            ro = sqrt(xo**2+yo**2)
            if(ro.lt.rstar) then
c	entire cell passes through star
              astar = one
            elseif((abs(x).lt.dx).and.(abs(y).lt.dx)) then
c	on one of inner 4 cells
              if(rstar.lt.dx) then
                astar = (pi/4.)*rstar**2/dx**2
              else
                astar = one-(sqrt2-rstar/dx)**2
              endif
            else
c	calculate fraction of cell passing through star
              ri = sqrt(x**2+y**2)-dx2
              ro = ri+dx
              if(ri.gt.rstar) then
                astar = zero
              elseif(ro.lt.rstar) then
                astar = one
              else
                astar = (rstar-ri)/dx
              endif
            endif
            if(astar.gt.zero) then
c	part of cell passed through star
              do iu = 1,nu
                fwi(iu) = fwi(iu)*(one-astar)
              enddo
              if((it.lt.1).or.(it.gt.nt)) return
              do iu = 1,nu
                ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi(iu)
              enddo
              fn(it,ir,ig) = fn(it,ir,ig)+gwt
              return
            endif
          else
            astar = zero
          endif
        elseif(r.lt.rstar) then
c	still in star off of midplane
          do iu = 1,nu
            fwi(iu) = zero
          enddo
          return
        endif

c	not on star
c	radiative transfer in cell
c	optical depth reduction factor is built into opa and sca
      if(opa(iuv,it,ir,ig)*ds.gt.tiny) then
        do iu = 1,nu
          fwi0 = fwi(iu)
          opac = opa(iu,it,ir,ig)
          scat = sca(iu,it,ir,ig)
          if(iy.le.ly) then
            sfun = swd(iu,it,ir,ig)
          else
            sfun = zero
          endif
          if((ig.gt.jg).and.(ig.le.ng)) then
            opac = dirt*opac+dir*opa(iu,it,jr,ig)+dit*opa(iu,jt,ir,ig)
            scat = dirt*scat+dir*sca(iu,it,jr,ig)+dit*sca(iu,jt,ir,ig)
            if(iy.le.ly)
     &        sfun = dirt*sfun+dir*swd(iu,it,jr,ig)+dit*swd(iu,jt,ir,ig)
          endif
          if(sfun.lt.zilch) sfun = zero
          opacs = opac+scat
          taud = opacs*ds
          if(taud.lt.small) then
            fwi1 = fwi0+(sfun-fwi0)*taud
            fwi2 = (fwi0+fwi1)/two
          else
c            expd = sexp(-taud)
            arg = hun*taud
            if(arg.lt.big) then
              iarg = int(arg)
              darg = (arg-iarg)/hun
              expd = (one-darg+darg**2/too)*expi(iarg+1)
              fwi1 = fwi0*expd+sfun*(one-expd)
            else
              fwi1 = sfun
            endif
c	Yes.  sfun does belong here.
            fwi2 = (fwi0-fwi1)/taud+sfun
          endif
          if(fwi1.lt.zilch) then
            fwi1 = zero
            fwi2 = fwi0/two
          endif
          fwi(iu) = fwi1
          ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi2
c	try using incoming radiation field?
c          ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi0
        enddo  !  iu
      else
        do iu = 1,nu
          ff(iu,it,ir,ig) = ff(iu,it,ir,ig)+gwt*fwi(iu)
        enddo
      endif

        fn(it,ir,ig) = fn(it,ir,ig)+gwt

        return
      end



      subroutine gastemp(it,ir,ig,iamax,tsubli,tgi,tdi)
c	calculate equilibrium gas and dust temperatures
c	holding pressure constant

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,epspe,agas,sigma,
     &          colin,colup,xco,xwat
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)
        common /taus/ taucont,tauline,tauterm,iuv,iusil,iufl,iwterm
        dimension tda(ma)
        logical doprint
        logical badnum

        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (fifth = 0.2d+00)
        parameter (small = 1.0d-03)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (five = 5.0d+00)
        parameter (ten = 1.0d+01)
        parameter (pi = 3.14159265d+00)
        parameter (bk = 1.3807d-16)
        parameter (sigsb = 5.670d-05)
        parameter (xmg = 4.0d-24)
        tpi = two*pi
        fpi = four*pi
        elog = exp(one)

        doprint = (myrank.eq.mysize/2)
        if(badnum(tgi).or.(tgi.gt.tstar)) tgi = tstar
        tgi = max(ten,tgi)
        dens0 = dg(it,ir,ig)
        pres = dens0*tgi
        dens = pres/tgi
        if(ig.eq.ng9) then
          dr = dr9
        else
          dr = dr0*two**ig
        endif
        r = dr*(ir+half)
        if((r.le.rstar).or.(dens.eq.zero)) goto 90
        if(nit.eq.-1) then
c	initialize opa and sca
          dd0 = dd(it,ir,ig)
          do iu = 1,nu
            opac = zero
            scat = zero
            do ia = 1,na
c	include settling factor in popa
              popa = condens(ia,it,ir,ig)*popd(ia)*sexp(-dd0*ad(ia))
     &               *pi*ad(ia)**2
              opac = opac+popa*phab(iu,ia)
              scat = scat+popa*phsc(iu,ia)
            enddo
            opa(iu,it,ir,ig) = opac*dens0
            sca(iu,it,ir,ig) = scat*dens0
          enddo
          dtauv = (opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
          taumx = quart*(taucont+half)
          if(dtauv.gt.taumx) then
            taufac = taumx/dtauv
            do iu = 1,nu
              opa(iu,it,ir,ig) = taufac*opa(iu,it,ir,ig)
              sca(iu,it,ir,ig) = taufac*sca(iu,it,ir,ig)
            enddo
          else
            taufac = one
          endif
          opl(ml,it,ir,ig) = taufac
        else
c	will divide taufac out of dust heating and cooling
          taufac = opl(ml,it,ir,ig)
        endif  !  first C&G iteration
        if(doprint.and.(taufac*(one+small-taufac).le.zero))
     &    print '("taufac",3i6,e10.2)', it,ir,ig,taufac
        densfac = dens*taufac
        if(densfac.le.zero) goto 90
c	heating and cooling powers/gas particle
        hdust = zero
        hpe = zero
        if(opa(iuv,it,ir,ig).gt.zero) then
          npe = nu-4
          do iu = nu,npe,-1
            opadu = opa(iu,it,ir,ig)*dud(iu)
            hdust = hdust+opadu*fpi*gg(iu,it,ir,ig)
          enddo
c	epspe should depend on dens and maybe ad
          hpe = epspe*hdust
          do iu = npe-1,1,-1
            opadu = opa(iu,it,ir,ig)*dud(iu)
            hdust = hdust+opadu*fpi*gg(iu,it,ir,ig)
          enddo
          hdust = hdust+opa(1,it,ir,ig)*ud(1)*gg(1,it,ir,ig)/four
        endif
c	divide out densfac to calculate heating and cooling per gas particle
c	note: densfac must use the dens used to calculate opa and sca
        hdust = hdust/densfac
        hpe = hpe/densfac
        hxr = (pxr/(fpi*r**2))*sexp(-axr*colin)*axr
        if(sigma.gt.zero) then
          haccr = gmmdot/(fpi*r**3*max(sigma,1.0d+18))
        else
          haccr = zero
        endif
        hgas = hpe+hxr+haccr
        htot = hgas+hdust
        popsum = zero
        dd0 = dd(it,ir,ig)
        do ia = 1,na
          popa = condens(ia,it,ir,ig)*popd(ia)*sexp(-dd0*ad(ia))
     &           *pi*ad(ia)**2
          popsum = popsum+popa
        enddo
c	molecular cooling function from Hollenbach & McKee (1979)
c	also calculates agas = dlnC/dlnT
c        cgas = cooling(dens,tgi,it,ir,ig)
c        if(popsum.eq.zero) then
c          if((hgas.gt.zero).and.(cgas.gt.zero).and.(agas.gt.one)) then
c            tgi = tgi*(hgas/cgas)**(one/agas)
c            dens = pres/tgi
c            cgas = cooling(dens,tgi,it,ir,ig)
c            tgi = tgi*(hgas/cgas)**(one/agas)
c            tgi = min(tstar,tgi)
c          else
c            tgi = 0.84*tstar*(rstar/r)**0.4
c          endif
c          tdi = tsubli
c          goto 90
c        endif
c	equilibrium temperature assuming tg = td
c        tgd = tgi
c        dens = pres/tgd
c        iit = 0
c  100   iit = iit+1
c          cgas = cooling(dens,tgd,it,ir,ig)
c          cdust = zero
c          do iu = 1,nu
c            opadu = opa(iu,it,ir,ig)*dud(iu)
c            cdust = cdust+opadu*fpi*bw(ud(iu),tgd)
c          enddo
c          cdust = cdust+opa(1,it,ir,ig)*ud(1)*fpi*bw(ud(1),tgd)/four
c          cdust = cdust/densfac
c          ctot = cgas+cdust
c          if(badnum(ctot).or.(ctot.le.zero)) then
c            if(doprint) then
c              print '("ctot",3i6,1p,5e10.2/3e10.2)',
c     &        it,ir,ig,ctot,cgas,cdust,dens,tgd,
c     &        opa(1,it,ir,ig),bw(ud(1),tgd),densfac
c            else
c              call sleep(1)
c            endif
c            stop
c          endif
c          atot = (agas*cgas+five*cdust)/ctot
cc	should this ever happen?
c          if(atot.le.zero) then
c            if(doprint) print '(e10.3,3i6,4e10.2)',
c     &        atot,it,ir,ig,agas,cgas,five,cdust
c            if(tgi.gt.tsubli/two) then
c              tgi  = tstar
c              tdi = tsubli
c            else
c              tdi = tgi
c            endif 
c            goto 90
c          elseif(atot.lt.half) then
c            atot = half
c          endif
c          hcrat = max(half,min(two,htot/ctot))
c          tgd = tgd*hcrat**(one/atot)
c          if(badnum(tgd).or.(tgd.le.zero)) then
c            print '("tgd",4i4,1p,6e10.2)',
c     &        it,ir,ig,iit,hdust,cdust,cgas,hcrat,atot,tgd
c            stop
c          endif
c          dens = pres/tgd
c	if(doprint.and.(it.eq.67).and.(ir.eq.67)) print '(i4,i3,1p,8e9.2)',
c     &	ig,iit,dens,tgi,hgas,cgas,agas,hdust,cdust,tgd
c        if((abs(hcrat-one).gt.small).and.(iit.lt.16).and.(tgd.lt.tstar))
c     &    goto 100
c        tgd = max(ten,min(tstar,tgd))
c        tdi = min(tgd,tsubli)
c        tgi = tgd
        dens = pres/tgi
        cgas = cooling(dens,tgi,it,ir,ig)
        acgd = four*dens*popsum*sqrt(bk*tgi/xmg)*acom*2.5*bk
c        if(tgd.ge.tstar) then
c          tgi = tstar
c        elseif(cdust.le.zero) then
c          goto 90
c        elseif((hgas.le.zero).and.(cgas.le.zero)) then
c          continue
c        elseif(acgd.le.zero) then
        if(acgd.le.zero) then
          if((hgas.gt.zero).and.(cgas.gt.zero)) then
            iit = 0
  112       iit = iit+1
              cgas = cooling(dens,tgi,it,ir,ig)
              tgi = tgi*(hgas/cgas)**(one/max(agas,half))
            if((iit.lt.16).and.(abs(hgas/cgas-one).gt.small)) goto 112
            tgi = min(tstar,max(ten,tgi))
          endif
          if(hdust.gt.zero) then
            iit = 0
  114       iit = iit+1
              cdust = zero
              do iu = 1,nu
                opadu = opa(iu,it,ir,ig)*dud(iu)
                cdust = cdust+opadu*fpi*bw(ud(iu),tdi)
              enddo
              cdust = cdust+opa(1,it,ir,ig)*ud(1)*fpi*bw(ud(1),tdi)/four
              cdust = cdust/densfac
              tdi = tdi*(hdust/cdust)**fifth
            if((iit.lt.16).and.(abs(hdust/cdust-one).gt.small)) goto 114
          else
            tdi = tgi
          endif
        else
c	calculate tg and td including gas-dust heat flow
          iit = 0
  110     iit = iit+1
            cdust = zero
            do iu = 1,nu
              opadu = opa(iu,it,ir,ig)*dud(iu)
              cdust = cdust+opadu*fpi*bw(ud(iu),tdi)
            enddo
            cdust = cdust+opa(1,it,ir,ig)*ud(1)*fpi*bw(ud(1),tdi)/four
            cdust = cdust/densfac
            cgas = cooling(dens,tgi,it,ir,ig)
            dtd = (hdust+acgd*(tgi-tdi)-cdust)/(five*cdust/tdi+acgd)
            agd = max((agas*cgas/tgi+acgd),half*cgas/tgi)
            dtg = (hgas+acgd*(tdi-tgi)-cgas)/agd
            dtmx = tdi/four
            tdi = tdi+min(dtmx,max(-dtmx,dtd))/two
            dtmx = tgi/four
            tgi = tgi+min(dtmx,max(-dtmx,dtg))/two
            tgi = min(tstar,max(ten,tgi))
            dens = pres/tgi
            if(badnum(tdi).or.badnum(tgi).or.(tgi.le.zero)) then
              print '("tgi",4i4,1p,6e10.2/7e10.2)',
     &        it,ir,ig,iit,hcdust,hcgas,dcdust,dcgas,acgd,dc2,
     &        dens,popsum,dtd,dtg,tdi,tgi
              goto 190
            endif
	if((ig.eq.ng9).and.(it.eq.nt).and.((ir-256)*(ir-320)*(ir-384).eq.0))
     &  print '(2i4,1p,6e10.2)', ir,iit,tgi,tdi,pres,dens,hgas,hdust
          if((abs(dtg/tgi).gt.small).and.(iit.lt.16)) goto 110
  190   endif
c       tg(it,ir,ig) = tgi  !  in main
c	dg(it,ir,ig) = dens  !  in hydro equil calc
c	calculate individual grain temperatures
c	using heat flows per grain surface area
        acgd = sqrt(bk*tgi/xmg)*acom*2.5*bk*dens
        iamax = 0
        do ia = 1,na
          if(condens(ia,it,ir,ig).gt.zero) then
            tdia = tdi
            hdia = zero
            do iu = 1,nu
              hdia = hdia+phab(iu,ia)*pi*gg(iu,it,ir,ig)*dud(iu)
            enddo
            iit = 0
  130       iit = iit+1
              cdia = zero
              do iu = 1,nu
                cdia = cdia+phab(iu,ia)*pi*bw(ud(iu),tdia)*dud(iu)
              enddo
              dcdia = five*cdia/tdia+acgd
              dtdia = (hdia-cdia+acgd*(tgi-tdia))/dcdia
              dtmx = tdia/four
              tdia = tdia+min(dtmx,max(-dtmx,dtdia))
              if(badnum(tdia).or.(tdia.le.zero)) then
                print '("tdia",5i4,1p,4e10.2/7e10.2)',
     &            it,ir,ig,ia,iit,hdia,cdia,dcdia,tdia,
     &              tgi,tdi,dens,acgd,hgas,cgas
                stop
              endif
            if((abs(dtdia/tdia).gt.small).and.(iit.lt.16)) goto 130
            tda(ia) = tdia
c	condens is now set in main
c	but should I let small grains sublimate here?
            if(tdia.gt.tsubli) then
              iamax = ia
            endif
          else
            tda(ia) = tdi
          endif
        enddo

        dd0 = dd(it,ir,ig)
        do iu = 1,nu
          opac = zero
          scat = zero
          sfun = zero
          do ia = 1,na
            condi = condens(ia,it,ir,ig)*sexp(-dd0*ad(ia))
            if(condi.gt.zero) then
              popa = condi*popd(ia)*pi*ad(ia)**2
              opac = opac+popa*phab(iu,ia)
              scat = scat+popa*phsc(iu,ia)
              sfun = sfun+popa*phab(iu,ia)*bw(ud(iu),tda(ia))
            endif
          enddo
          opa(iu,it,ir,ig) = opac*dens
          sca(iu,it,ir,ig) = scat*dens
          if((opac+scat).gt.zero) then
            swd(iu,it,ir,ig) =
     &        (sfun+scat*gg(iu,it,ir,ig))/(opac+scat)
          else
            swd(iu,it,ir,ig) = gg(iu,it,ir,ig)
          endif
        enddo
        dtauv = (opa(iuv,it,ir,ig)+sca(iuv,it,ir,ig))*dr
        taumx = quart*(taucont+fifth)
        if(dtauv.gt.taumx) then
          taufac = taumx/dtauv
        else
          taufac = one
        endif
        opl(ml,it,ir,ig) = taufac
        do iu = 1,nu
          if(taufac.lt.one) then
            opa(iu,it,ir,ig) = taufac*opa(iu,it,ir,ig)
            sca(iu,it,ir,ig) = taufac*sca(iu,it,ir,ig)
          endif
        enddo

        return

   90   tdi = min(tstar,tdi)
        tgi = min(tstar,max(ten,tgi))
        do iu = 1,nu
          w = ud(iu)
          opa(iu,it,ir,ig) = zero
          sca(iu,it,ir,ig) = zero
          swd(iu,it,ir,ig) = bw(w,tdi)
          opl(ml,it,ir,ig) = one
        enddo
        do ia = 1,na
          condens(ia,it,ir,ig) = zero
        enddo
        iamax = na+1
        return

      end



      function cooling(dg0,tg0,it,ir,ig)
c	gas cooling by H2O and CO rot lines from Hollenbach & McKee (1979)
c	per gas particle, T dependance at constant pressure

        implicit real*8 (a-h,o-z)
        parameter (mt = 192)
        parameter (mr = 128)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (mu = 64)
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,epspe,agas,sigma,
     &          colin,colup,xco,xwat
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        parameter (ten = 1.0d+01)
        parameter (pi = 3.14159265d+00)
        parameter (bk = 1.3807d-16)
        parameter (hck = 1.438d+00)
        tpi = two*pi
        fpi = four*pi
        elog = exp(one)
        deltat = 1.0d-03*tg0

        dvturb = dvturb0
        pres = dg0*tg0
        if(pres*(xco+xwat).le.zero) then
          cooling = zero
          agas = one
          return
        endif
        if(colin*colup.gt.zero) then
          colg = one/(one/colin+one/colup)
        else
          colg = zero
        endif
        ds = colg/dg0

c	water
        if(xwat.gt.zero) then
          amol = 6.0d-03
          tfreez = 1.0d+02
          trot0 = 2.3d+01
          dcrit = 7.3d+08
          iumol = 16.*slog10(sqrt(trot0*tg0)/hck)-16.
          tgas = tg0
          dens = pres/tgas
          xmol = xwat*sexp(-(tfreez/tgas)**10)
          colmol = xwat*colg
          taudmol = opa(iumol,it,ir,ig)*ds
          eta = sqrt(one+0.06d+00*tgas)
          coltau = 1.18d+07*dvturb*tgas*trot0**2/amol
          taut = four*colmol/(eta*coltau)
          ctau = taut*sqrt(tpi*slog(2.13d+00+(taut/elog)**2))
          dcritd = dcrit*sqrt(tgas)/dens
          ym = slog(one+ctau/(one+ten*dcritd))
          cwat = xmol*(two*bk*amol/trot0)*(two+ym+0.6d+00*ym**2)
     &      *sexp(-taudmol)/(one+ctau+dcritd+1.5d+00*sqrt(dcritd))
c	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          xmol = xwat*sexp(-(tfreez/tgas)**10)
          eta = max(one,0.24d+00*sqrt(tgas))
          coltau = 1.18d+07*dvturb*tgas*trot0**2/amol
          taut = four*colmol/(eta*coltau)
          ctau = taut*sqrt(tpi*slog(2.13d+00+(taut/elog)**2))
          dcritd = dcrit*sqrt(tgas)/dens
          ym = slog(one+ctau/(one+ten*dcritd))
          cwat1 = xmol*(two*bk*amol/trot0)*(two+ym+0.6d+00*ym**2)
     &      *sexp(-taudmol)/(one+ctau+dcritd+1.5d+00*sqrt(dcritd))
        else
          cwat = zero
          cwat1 = zero
        endif

c	CO
        if(xco.gt.zero) then
          amol = 9.7d-08
          tfreez = 3.0d+01
          trot0 = 2.76d+00
          dcrit = 1.9d+04
          iumol = 16.*slog10(sqrt(trot0*tg0)/hck)-16.
          tgas = tg0
          dens = pres/tgas
          xmol = xco*sexp(-(tfreez/tgas)**10)
          colmol = xco*colg
          taudmol = opa(iumol,it,ir,ig)*ds
          eta = one
          coltau = 1.18d+07*dvturb*tgas*trot0**2/amol
          taut = four*colmol/(eta*coltau)
          ctau = taut*sqrt(tpi*slog(2.13d+00+(taut/elog)**2))
          dcritd = dcrit*sqrt(tgas)/dens
          ym = slog(one+ctau/(one+ten*dcritd))
          cco = xmol*(two*bk*amol/trot0)*(two+ym+0.6d+00*ym**2)
     &      *sexp(-taudmol)/(one+ctau+dcritd+1.5d+00*sqrt(dcritd))
c	repeat with 1.001*tgas for T dependance
          tgas = tg0+deltat
          dens = pres/tgas
          xmol = xco*sexp(-(tfreez/tgas)**10)
          coltau = 1.18d+07*dvturb*tgas*trot0**2/amol
          taut = four*colmol/(eta*coltau)
          ctau = taut*sqrt(tpi*slog(2.13d+00+(taut/elog)**2))
          dcritd = dcrit*sqrt(tgas)/dens
          ym = slog(one+ctau/(one+ten*dcritd))
          cco1 = xmol*(two*bk*amol/trot0)*(two+ym+0.6d+00*ym**2)
     &      *sexp(-taudmol)/(one+ctau+dcritd+1.5d+00*sqrt(dcritd))
        else
          cco = zero
          cco1 = zero
        endif

c	H2 wild guess
        trot0 = 3.5d+02
        iumol = 16.*slog10(sqrt(trot0*tg0)/hck)-16.
        taudmol = opa(iumol,it,ir,ig)*ds
        chh0 = 1.0d-35
        tgas = tg0
        chh = chh0*tgas**3*sexp(-taudmol)
        tgas = tg0+deltat
        chh1 = chh0*tgas**3*sexp(-taudmol)

        cooling = cwat+cco+chh
        cool1 = cwat1+cco1+chh1
        if(cooling*deltat.gt.zero) then
          agas = (tg0/cooling)*(cool1-cooling)/deltat
        else
          agas = one
        endif

        return
      end



      subroutine ggavg(nul,ntsm,nrsm)

        implicit real*8 (a-h,o-z)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mg = 20)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)
        dimension ggsm(mu,mt),gnsm(mt)
        logical doprint,even
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)

        doprint = (myrank.eq.mysize/2)
        nr2 = nr/2

        if(ntsm.gt.1) then
          nsm2 = (ntsm-1)/2
          even = (nsm2.ne.(ntsm/2))
c	include angular smoothing in grid 9, but theta_sm will be smaller
          do igr = 0,ng
c	do only 2 point smooth in grids 1-8
            if(igr.gt.0) then
              nsm2 = 0
              even = .true.
            endif
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            xlr = lr
            do ir = 1,lr
              do it = 1,lt
                jt1 = it-nsm2
                jt2 = it+nsm2
                gnsm(it) = zero
                do iu = 1,nul
                  ggsm(iu,it) = zero
                enddo
                do jt = jt1,jt2
                  if(jt.lt.1) then
                    kt = 1
                  elseif(jt.gt.nt) then
                    kt = 2*nt-jt+1
                  elseif(jt.gt.lt) then
                    kt = lt
                  else
                    kt = jt
                  endif
                  gnsm(it) = gnsm(it)+gn(kt,ir,ig)
                  do iu = 1,nul
                    ggsm(iu,it) = ggsm(iu,it)+gg(iu,kt,ir,ig)
                  enddo
                enddo
                if(even) then
                  jt1 = jt1-1
                  jt2 = jt2+1
                  if(jt1.lt.1) then
                    jt1 = 1
                  elseif(jt2.gt.nt) then
                    jt2 = 2*nt-jt2+1
                  elseif(jt2.gt.lt) then
                    jt2 = lt
                  endif
                  gnsm(it) = gnsm(it)+(gn(jt1,ir,ig)+gn(jt2,ir,ig))/two
                  do iu = 1,nul
                    ggsm(iu,it) = ggsm(iu,it)
     &                +(gg(iu,jt1,ir,ig)+gg(iu,jt2,ir,ig))/two
                  enddo
                endif
              enddo  !  it
              rrat = (xlr/ir)**2
              do it = 1,lt
                gn(it,ir,ig) = rrat*gnsm(it)
                do iu = 1,nul
                  gg(iu,it,ir,ig) = rrat*ggsm(iu,it)
                enddo
              enddo  !  it
            enddo  !  ir
          enddo  !  igr
        endif

c	smooth mean intensities radially 2 points
        if(nrsm.gt.1) then
          xn9 = two/(one+two*nn9)
          kr2 = nr/2+1
          do igr = 0,ng
            call startgr(igr,ig,kr,lr,lt,dr,dcosth)
            kr = max(kr,2)
            do it = 1,lt
              it1 = nt-(nt-it)/nn9
              do iu = 1,nul
                ggj = gg(iu,it,kr-1,ig)
                ggk = gg(iu,it,kr,ig)
                do ir = kr,lr-1
                  ggi = ggj
                  ggj = ggk
                  ggk = gg(iu,it,ir+1,ig)
                  gg(iu,it,ir,ig) = (ggi+two*ggj+ggk)/four
                enddo
                if(ig.lt.ng) then
                  gg(iu,it,lr,ig) = ((ggi+ggj)/two+two*ggk
     &              +gg(iu,it,kr2,ig+1))/four
                elseif(ig.eq.ng9) then
                  ggi = gg(iu,it1,kr2,2)
                  gg(iu,it,lr,ig) = (ggj+(3.-xn9)*ggk+xn9*ggi)/four
                endif
              enddo
              gnj = gn(it,kr-1,ig)
              gnk = gn(it,kr,ig)
              do ir = kr,lr-1
                gni = gnj
                gnj = gnk
                gnk = gn(it,ir+1,ig)
                gn(it,ir,ig) = (gni+two*gnj+gnk)/four
              enddo
              if(ig.lt.ng) then
                gn(it,lr,ig) = ((gni+gnj)/two+two*gnk
     &            +gn(it,kr2,ig+1))/four
              elseif(ig.eq.ng9) then
                gni = gn(it1,kr2,2)
                gn(it,lr,ig) = (gnj+(3.-xn9)*gnk+xn9*gni)/four
              endif
            enddo  !  it
          enddo  !  igr
        endif

c	divide gg by gn
        gnmin = 1.0d+12
        gnmax = zero
        gnsum = zero
        sumi = zero
        nzero = 0
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          ls = nt-lt+1
          do ir = kr,lr
            r = dr*(ir-half)
            do it = 1,lt
              gni = gn(it,ir,ig)
              if(gni.eq.zero) then
                if((r.gt.rstar).and.(opl(1,it,ir,ig).gt.zero)) then
                  if(doprint.and.(nzero.eq.0))
     &              print'("Warning: no rays pass through ig,ir,it =",
     &              3i6)', ig,ir,it
                  nzero = nzero+1
                endif
              elseif((ir.gt.lr/2).and.(it.gt.(nt/2))) then
                if(gni.lt.gnmin) then
                  gnmin = gni
                  itmin = it
                  irmin = ir
                  igmin = ig
                elseif(gni.gt.gnmax) then
                  gnmax = gni
                  itmax = it
                  irmax = ir
                  igmax = ig
                endif
                gnsum = gnsum+gni
                sumi = sumi+one
              endif
              if(gni.gt.zero) then
                do iu = 1,nul
                  gg(iu,it,ir,ig) = gg(iu,it,ir,ig)/gni
                enddo
              elseif(ir.gt.kr) then
                do iu = 1,nul
                  gg(iu,it,ir,ig) = gg(iu,it,ir-1,ig)
                enddo
              else
                do iu = 1,nul
                  gg(iu,it,ir,ig) = zero
                enddo
              endif
            enddo  !  it
          enddo  !  ir
        enddo  !  igr
        gnmean = gnsum/sumi
        if(gnmin.gt.gnmax) gnmin = zero
        if(doprint.and.(nal.eq.na)) then
          print '(" min,max,mean rays per cell:",1p,3f12.2)',
     &      gnmin,gnmax,gnmean
          print '(6i8)', itmin,irmin,igmin,itmax,irmax,igmax
          if(nzero.gt.0)
     &    print '(" no rays pass through",i8," cells of",i8)',
     &      nzero,nt*nr2*(ng-1)+nt*nr9
        endif

        return
      end



      subroutine grid0(fwk,ly)

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwk(mw)
        dimension fwi(mw)
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /angls/ dcost,dcos9,rot,irot
        common /taus/ taucont,tauline,tauterm,iuv,iusil,iufl,iwterm
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

c	first grid to calculate
        ig = ng
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

c	zero spectrum summed over lines of sight
        do iw = 1,lw
          fwk(iw) = zero
        enddo
        do ix = 1,nx
          do iy = 1,ly
            taucont = zero
            tauline = zero
            tauterm = zero
c	in middle 1/2 x 1/2 of grid?
            middle = (ix.gt.nx4).and.(ix.le.nx3)
     &        .and.(iy.gt.nx4).and.(iy.le.nx3)
            if(middle) then
c	use finer grid for grid section above and below middle
              jx2 = 2*ix-nx2
              jx1 = jx2-1
              jy2 = 2*iy-nx2
              jy1 = jy2-1
              do jx = jx1,jx2
                do jy = jy1,jy2
                  taucont = zero
                  tauline = zero
                  tauterm = zero
c	zero spectrum for line of sight
                  do iw = 1,lw
                    fwi(iw) = zero
                  enddo
                  iz = 0
                  wdopp = zero
c	step along line of sight
  14              iz = iz+1
                  if(iz.eq.nz5) then
                    if(ig.eq.12) then
                      call gridb(fwi,jx,jy,wdopp)
                    elseif(ig.eq.11) then
                      call grida(fwi,jx,jy,wdopp)
                    elseif(ig.eq.10) then
                      call grid9(fwi,jx,jy,wdopp)
                    elseif(ig.eq.9) then
                      call grid8(fwi,jx,jy,wdopp)
                    elseif(ig.eq.8) then
                      call grid7(fwi,jx,jy,wdopp)
                    elseif(ig.eq.7) then
                      call grid6(fwi,jx,jy,wdopp)
                    elseif(ig.eq.6) then
                      call grid5(fwi,jx,jy,wdopp)
                    elseif(ig.eq.5) then
                      call grid4(fwi,jx,jy,wdopp)
                    elseif(ig.eq.4) then
                      call grid3(fwi,jx,jy,wdopp)
                    elseif(ig.eq.3) then
                      call grid2(fwi,jx,jy,wdopp)
                    else
                      print '(" ng =",i4,"?")', ig
                      stop
                    endif
                    iz = iz+nz/2
                  endif
c	call radiative transfer routine for step
                  call step(fwi,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
                  if(iz.lt.nz) go to 14
c	store compressed spectrum in data cube
                  if(ig.le.ngc) call storecube(fwi,jg,jx,jy)
c	add outgoing spectrum to sum from grid
                  do iw = 1,lw
                    fwk(iw) = fwk(iw)+fwi(iw)/4.
                  enddo
c	end loops over 4 lines of sight
                enddo
              enddo
            else
c	not in middle.  follow ray along ix,iy
              do iw = 1,lw
                fwi(iw) = zero
              enddo
              wdopp = zero
c	step along line of sight
              do iz = 1,nz
c	call radiative transfer routine for step
                call step(fwi,wdopp,ix,iy,iz,ig)
              enddo
c	store compressed spectrum in data cube
              if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
c	add outgoing spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwi(iw)
              enddo
            endif
c	end loops over lines of sight in ix,iy grid
          enddo
        enddo

        return
      end


      subroutine gridb(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 11
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grida(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grida(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 10
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid9(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid9(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 9
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid8(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid8(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 8
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid7(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid7(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 7
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid6(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid6(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 6
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid5(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid5(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 5
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid4(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid4(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 4
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid3(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
c	call radiative transfer routine for step
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid3(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 3
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid2(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
c	call radiative transfer routine for step
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid2(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = 2
        jg = ig-1
        nx2 = nx/2
        nx4 = nx/4
        nx3 = 3*nx4
        nz5 = nz/4+1

        middle = (ix.gt.nx4).and.(ix.le.nx3)
     &    .and.(iy.gt.nx4).and.(iy.le.nx3)
        if(middle) then
          do iw = 1,lw
            fwk(iw) = zero
          enddo
c	branch out by 2x2
          jx2 = 2*ix-nx2
          jx1 = jx2-1
          jy2 = 2*iy-nx2
          jy1 = jy2-1
          do jx = jx1,jx2
            do jy = jy1,jy2
c	copy spectrum from outer grid
              do iw = 1,lw
                fwj(iw) = fwi(iw)
              enddo
              wdopp = wdopp0
c	step along line of sight
              iz = 0
  24          iz = iz+1
              if(iz.eq.nz5) then
                call grid1(fwj,jx,jy,wdopp)
                iz = iz+nz/2
              endif
c	call radiative transfer routine for step
              call step(fwj,wdopp,jx,jy,iz,-ig)
c	end loop along line of sight
              if(iz.lt.nz) go to 24
c	store compressed spectrum in data cube
              if(ig.le.ngc) call storecube(fwj,jg,jx,jy)
c	add spectrum to sum from grid
              do iw = 1,lw
                fwk(iw) = fwk(iw)+fwj(iw)/4.
              enddo
c	end loops over 4 lines of sight
            enddo
          enddo
          do iw = 1,lw
            fwi(iw) = fwk(iw)
          enddo
        else
c	not in the middle, don''t branch out
          wdopp = wdopp0
c	step along line of sight
          do iz = 1,nz
c	call radiative transfer routine for step
            call step(fwi,wdopp,ix,iy,iz,ig)
          enddo
c	store compressed spectrum in data cube
          if(ig.eq.ngc) call storecube(fwi,ig,ix,iy)
        endif
        wdopp0 = wdopp

        return
      end


      subroutine grid1(fwi,ix,iy,wdopp0)

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mu = 64)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        dimension fwi(mw)
        dimension fwj(mw),fwk(mw)
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /angls/ dcost,dcos9,rot,irot
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        ig = ng9
c	middle is not a special case for grid 1
c	store expanded dust grid to cube(,,,2)

        do iw = 1,lw
          fwk(iw) = zero
        enddo
c	branch out by nn9*nn9
        jx1 = nn9*(ix-1)+1
        jx2 = nn9*ix
        jy1 = nn9*(iy-1)+1
        jy2 = nn9*iy
        do jx = jx1,jx2
          do jy = jy1,jy2
c	copy spectrum from outer grid
            do iw = 1,lw
              fwj(iw) = fwi(iw)
            enddo
            wdopp = wdopp0
c	step along line of sight
            do iz = 1,nz9
              call step(fwj,wdopp,jx,jy,iz,ig)
            enddo
c	store compressed dust spectrum in data cube
            if(moltype.le.0) call storecube(fwj,ig,jx,jy)
c	add spectrum to sum from grid
            do iw = 1,lw
              fwk(iw) = fwk(iw)+fwj(iw)
            enddo
c	end loops over lines of sight on grid
          enddo
        enddo
c	put sum from grid into fwi
        xn9 = nn9**2
        do iw = 1,lw
          fwi(iw) = fwk(iw)/xn9
        enddo
        wdopp0 = wdopp

        return
        end



      subroutine storecube(fwi,ig,ix,iy)
c	store compressed spectrum in data cube

        implicit real*8 (a-h,o-z)
        parameter (mw = 1024)
        parameter (ml = 2)
        parameter (mc = 256)
        dimension fwi(mw)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /taus/ taucont,tauline,tauterm,iuv,iusil,iufl,iwterm
        common /cube/ cube(mc,mc,mc),ngc
        logical middle

        if(ig.eq.ng9) then
          if(moltype.gt.0) return
          nc = mc
          if(nc.ge.nx9) then
            ncnx = nc/nx9
            kx4 = ncnx*ix
            kx3 = kx4-ncnx+1
            ky4 = ncnx*iy
            ky3 = ky4-ncnx+1
            if((kx3.lt.1).or.(kx4.gt.mc).or.(ky3.lt.1).or.(ky4.gt.mc)) then
              print '(" in storecube kx3,kx4,ky3,ky4 =",5i8)',
     &          kx3,kx4,ky3,ky4,mc
              return
            endif
            fwisil = fwi(iusil)
            do ky = ky3,ky4
              ky1 = nc+1-ky
              do kx = kx3,kx4
                cube(kx,ky,2) = fwisil
                cube(kx,ky1,2) = fwisil
              enddo
            enddo
          else
            ncnx = nx9/nc
            kx = (ix+1)/ncnx
            ky = (iy+1)/ncnx
            if((kx.lt.1).or.(kx.gt.nc).or.(ky.lt.1).or.(ky.gt.nc)) then
              print '(" in storecube kx3,kx4,ky3,ky4 =",5i8)',
     &          kx3,kx4,ky3,ky4,nc
              return
            endif
            fwisil = fwi(iusil)/ncnx
            ky1 = nc+1-ky
            if((mod(ix,ncnx).eq.1).and.(mod(iy,ncnx).eq.1)) then
              cube(kx,ky,2) = fwisil
              cube(kx,ky1,2) = fwisil
            else
              cube(kx,ky,2) = cube(kx,ky,2)+fwisil
              cube(kx,ky1,2) = cube(kx,ky1,2)+fwisil
            endif
          endif
        elseif(ig.gt.ngc) then
          return
        else
          nc = mc/2**(ngc-ig)
          ncnx = nc/nx
          middle = (ix.gt.nx/4).and.(ix.le.3*nx/4)
     &      .and.(iy.gt.nx/4).and.(iy.le.3*nx/4)
          if((middle.and.(ncnx.gt.1).and.(ig.gt.1))
     &      .or.(ncnx.lt.1)) return
          kx4 = ncnx*ix+(mc-nc)/2
          kx3 = kx4-ncnx+1
          ky4 = ncnx*iy+(mc-nc)/2
          ky3 = ky4-ncnx+1
          if((kx3.lt.1).or.(kx4.gt.mc).or.(ky3.lt.1).or.(ky4.gt.mc)) then
            print '(" in storecube kx3,kx4,ky3,ky4 =",5i8)',
     &        kx3,kx4,ky3,ky4,mc
            stop
          endif
          fwisil = fwi(iusil)
          do ky = ky3,ky4
            ky1 = mc+1-ky
            do kx = kx3,kx4
              cube(kx,ky,1) = fwisil
              cube(kx,ky1,1) = fwisil
            enddo
          enddo
        endif

        return
      end



      subroutine step(fwi,wdopp0,ix,iy,iz,igr)

        implicit real*8 (a-h,o-z)
        parameter (ml = 2)
        parameter (mw = 1024)

        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /angls/ dcost,dcos9,rot,irot
        dimension fwi(mw)

        call stepd(fwi,ix,iy,iz,igr)

        return
      end




      subroutine movie(nc,ic)

        implicit real*8 (a-h,o-z)
        parameter (mc = 256)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /angls/ dcost,dcos9,rot,irot
        common /cube/ cube(mc,mc,mc),ngc
        dimension frame(mc,mc)
        character*16 lockline

c	make and plot continuum frame
        cosrot = cos(rot)
        xmid = (nc+1)/two
        nc2 = nc/2
        sumcon = zero
        do iy = 1,nc
          do ix = 1,nc
            x = (ix-xmid)/cosrot
            y = iy-xmid
            scale = sqrt(x**2+y**2)
            scale = (x**2+y**2)
            if(scale.lt.one) scale = zero
            frame(ix,iy) = sqrt(max(zero,scale*cube(ix,iy,ic)))
            sumcon = sumcon+frame(ix,iy)
          enddo
        enddo
        if(sumcon.eq.zero) then
          print '(" continuum = 0")'
          return
        endif
   18   open(unit=8,file='plotlock',access='sequential',err=19)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 18
        endif
        open(unit=8,file='plotlock',access='sequential',err=19)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='images',access='direct',recl=2*nc)
        irec = 0
        do iy = 1,nc
          irec = irec+1
          write(unit=9,rec=irec) (frame(ix,iy),ix=1,nc)
        enddo
        nfr = 1

        close(unit=9)
        open(unit=8,file='plotlock',access='sequential',err=19)
        write(8,'("images")')
        write(8,'(3i8)') nc,nc,nfr
        close(unit=8)

   19   return
      end



      subroutine zerof(nul)
c	zero ff and fn

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ml = 2)
        parameter (mu = 64)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        parameter (zero = 0.0d+00)

        do ig = 1,ng
          do ir = 1,nr
            do it = 1,nt
              do il = 1,nul
                ff(il,it,ir,ig) = zero
              enddo
              fn(it,ir,ig) = zero
            enddo
          enddo
        enddo
        ig = ng9
          do ir = 1,nr9
            do it = 1,nt
              do il = 1,nul
                ff(il,it,ir,ig) = zero
              enddo
              fn(it,ir,ig) = zero
            enddo
          enddo

        return
      end


      subroutine zerog(nul)
c	zero gg and gn

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mu = 64)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /mean/ gg(mu,mt,mr,mg),gn(mt,mr,mg),gguv(mt,mr,mg),
     &          ff(mu,mt,mr,mg),fn(mt,mr,mg)
        parameter (zero = 0.0d+00)

        do ig = 1,ng
          do ir = 1,nr
            do it = 1,nt
              do il = 1,nul
                gg(il,it,ir,ig) = zero
              enddo
              gn(it,ir,ig) = zero
            enddo
          enddo
        enddo
        ig = ng9
          do ir = 1,nr9
            do it = 1,nt
              do il = 1,nul
                gg(il,it,ir,ig) = zero
              enddo
              gn(it,ir,ig) = zero
            enddo
          enddo

        return
      end



      subroutine copyarr(arr)
c	copy array between grids

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension arr(mt,mr,mg)
        logical badnum,doprint
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2
c	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        xn9 = nn9**2
        nn8 = nn9-1
        nt8 = nt-nt/nn9+1
        do ir = 1,nr
          jr2 = nn9*ir
          jr1 = jr2-nn8
          do it = 1,nt,nn9
            suma = zero
            do jr = jr1,jr2
              do jt = it,it+nn8
                suma = suma+arr(jt,jr,jg)
              enddo
            enddo
            if(badnum(suma)) then
              if(doprint) print '(e10.2," in copyarr at",3i6)',
     &          suma,it,jr1,jg
              suma = zero
            endif
            arr(it/nn9+nt8,ir,ig) = suma/xn9
          enddo
        enddo
c	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr2 = 2*ir
            jr1 = jr2-1
            do it = 1,nt
              arr(it,ir,ig) = (arr(it,jr1,jg)+arr(it,jr2,jg))/two
              if(badnum(arr(it,ir,ig))) then
                if(doprint) print '(e10.2," in copyarr at",3i6)',
     &            arr(it,ir,ig),it,jr1,jg
                arr(it,ir,ig) = zero
              endif
            enddo
          enddo
        enddo

        return
      end



      subroutine copytau(arr)
c	copy array between grids without averaging

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension arr(mt,mr,mg)
        logical badnum,doprint
        parameter (zero = 0.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2
c	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        xn9 = nn9
        nn8 = nn9-1
        nt8 = nt-nt/nn9+1
        do ir = 1,nr
          jr = nn9*ir
          do it = 1,nt,nn9
            suma = zero
            do jt = it,it+nn8
              suma = suma+arr(jt,jr,jg)
            enddo
            if(badnum(suma)) then
              if(doprint) print '(e10.2," in copytau at",3i6)',
     &          suma,it,jr,jg
              suma = zero
            endif
            arr(it/nn9+nt8,ir,ig) = suma/xn9
          enddo
        enddo
c	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr = 2*ir
            do it = 1,nt
              arr(it,ir,ig) = arr(it,jr,jg)
            enddo
          enddo
        enddo

        return
      end



      subroutine copywta(arr,wtarr)
c	copy between grids with weighting

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension arr(mt,mr,mg),wtarr(mt,mr,mg)
        logical badnum,doprint

        doprint = (myrank.eq.0)
        nr2 = nr/2
        nn8 = nn9-1
        nt8 = nt-nt/nn9+1
        xn9 = nn9**2
c	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        do ir = 1,nr
          jr2 = nn9*ir
          jr1 = jr2-nn8
          do it = 1,nt,nn9
            suma = zero
            sumwt = zero
            do jr = jr1,jr2
              do jt = it,it+nn8
                suma = suma+arr(jt,jr,jg)*wtarr(jt,jr,jg)
                sumwt = sumwt+wtarr(jt,jr,jg)
              enddo
              if(badnum(suma)) then
                if(doprint) print '(e10.2," in copywta at",3i6)',
     &            suma,it,jr1,jg
                suma = zero
              endif
            enddo
            if(sumwt.gt.zero) then
              arr(it/nn9+nt8,ir,ig) = suma/sumwt
            else
              suma = zero
              do jr = jr1,jr2
                do jt = it,it+nn8
                  suma = suma+arr(jt,jr,jg)
                enddo
                if(badnum(suma)) then
                  if(doprint) print '(e10.2," in copywta at",3i6)',
     &              suma,it,jr1,jg
                  suma = zero
                endif
              enddo
              arr(it/nn9+nt8,ir,ig) = suma/xn9
            endif
          enddo
        enddo
c	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr2 = 2*ir
            jr1 = jr2-1
            do it = 1,nt
              sumwt = wtarr(it,jr1,jg)+wtarr(it,jr2,jg)
              if(sumwt.gt.zero) then
                arr(it,ir,ig) = (arr(it,jr1,jg)*wtarr(it,jr1,jg)
     &            +arr(it,jr2,jg)*wtarr(it,jr2,jg))/sumwt
              else
                arr(it,ir,ig) = (arr(it,jr1,jg)+arr(it,jr2,jg))/two
              endif
              if(badnum(arr(it,ir,ig))) then
                if(doprint) print '(e10.2," in copywta at",3i6)',
     &            arr(it,ir,ig),it,jr1,jg
                arr(it,ir,ig) = zero
              endif
            enddo
          enddo
        enddo

        return
      end



      subroutine copyarrr(arrr,mul,nul)
c	copy 4-d array between grids

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        logical badnum,doprint
        dimension arrr(mul,mt,mr,mg)
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)

        doprint = (myrank.eq.0)
        nr2 = nr/2
        nn8 = nn9-1
        xn9 = nn9**2
        nt8 = nt-nt/nn9+1
c	copy from expanded grid onto grid 1
        ig = 1
        jg = ng9
        do ir = 1,nr
          jr2 = nn9*ir
          jr1 = jr2-nn8
          do it = 1,nt,nn9
            do il = 1,nul
              suma = zero
              do jr = jr1,jr2
                do jt = it,it+nn8
                  suma = suma+arrr(il,jt,jr,jg)
                enddo
              enddo
              if(badnum(suma)) then
                if(doprint) print '(e10.2," in copyarrr at",4i6)',
     &            suma,il,it,jr1,jg
                suma = zero
              endif
              arrr(il,it/nn9+nt8,ir,ig) = suma/xn9
            enddo
          enddo
        enddo
c	copy to inner half of next grid
        do ig = 2,ng
          jg = ig-1
          do ir = 1,nr2
            jr2 = 2*ir
            jr1 = jr2-1
            do it = 1,nt
              do il = 1,nul
                arrr(il,it,ir,ig)
     &            = (arrr(il,it,jr1,jg)+arrr(il,it,jr2,jg))/two
                if(badnum(arrr(il,it,ir,ig))) then
                  if(doprint) print '(e10.2," in copyarrr at",4i6)',
     &              arrr(il,it,ir,ig),il,it,jr1,jg
                  arrr(il,it,ir,ig) = zero
                endif
              enddo
            enddo
          enddo
        enddo

        return
      end



      subroutine plotcon(con,logplot,monfile)
c	make r-th contour plots

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (zero = 0.0d+00)
        parameter (two = 2.0d+00)
        parameter (four = 4.0d+00)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        dimension con(mt,mr,mg)
        logical logplot,monplot
        character*16 lockline,monfile
        dimension conmap(2*mt)

        monplot = (index(monfile,' ').gt.1)
        nr2 = nr/2

   86   open(unit=8,file='plotlock',access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 86
        endif
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='contours',access='direct',recl=2*nt)
        if(monplot) open(unit=10,file=monfile,access='direct',recl=2*nt)
 
c	plot map of grid ng9
        irec = 0
        ig = ng9
        if((nr9/2).le.nt) then
          nxm = nt
          nym = nr9/2
          do ir = nr9/2+1,nr9
            imap = 0
            do it = 1,nt
              imap = imap+1
              if(logplot) then
                if(con(it,ir,ig).gt.zero) then
                  conmap(imap) = slog10(con(it,ir,ig))
                else
                  conmap(imap) = zero
                endif
              else
                conmap(imap) = con(it,ir,ig)
              endif
            enddo  !  it
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=1,nt)
          enddo  !  ir
        elseif((nr9/2).le.(2*nt)) then
          nxm = 2*nt
          nym = nr9/2
          do ir = nr9/2+1,nr9
            imap = 0
            do it = 1,nt
              imap = imap+2
              if(logplot) then
                if(con(it,ir,ig).gt.zero) then
                  conmap(imap) = slog10(con(it,ir,ig))
                else
                  conmap(imap) = zero
                endif
              else
                conmap(imap) = con(it,ir,ig)
              endif
              if(imap.gt.2) then
                conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
              else
                conmap(1) = conmap(2)
              endif
            enddo  !  it
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=1,nt)
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
            if(monplot) then
              write(unit=10,rec=irec-1) (conmap(it),it=1,nt)
              write(unit=10,rec=irec) (conmap(it),it=nt+1,2*nt)
            endif
          enddo  !  ir
        else
          nxm = 2*nt
          nym = nr9/4
          do ir = nr9/2+1,nr9,2
            imap = 0
            do it = 1,nt
              imap = imap+2
              if(logplot) then
                prod = con(it,ir,ig)*con(it,ir+1,ig)
                if(prod.gt.zero) then
                  conmap(imap) = slog10(prod)/two
                else
                  conmap(imap) = zero
                endif
              else
                conmap(imap) = (con(it,ir,ig)+con(it,ir+1,ig))/two
              endif
              if(imap.gt.2) then
                conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
              else
                conmap(1) = conmap(2)
              endif
            enddo  !  it
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=1,nt)
            irec = irec+1
            write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
            if(monplot) then
              write(unit=10,rec=irec-1) (conmap(it),it=1,nt)
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
            endif
          enddo  !  ir
        endif

        close(unit=9)
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("contours")')
        write(8,'(3i8)') nxm,nym,1
        close(unit=8)
        call waitasec(.false.,4)

   88   open(unit=8,file='plotlock',access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 88
        endif
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='contours',access='direct',recl=2*nt)
 
c	plot map of grids 1-ng
        jrec = irec
        irec = 0
        if((ng*nr/2).le.(2*nt)) then
          nxm = 2*nt
          nym = ng*nr2
          do ig = 1,ng
            do ir = nr2+1,nr
              imap = 0
              do it = 1,nt
                imap = imap+2
                if(logplot) then
                  if(con(it,ir,ig).gt.zero) then
                    conmap(imap) = slog10(con(it,ir,ig))
                  else
                    conmap(imap) = zero
                  endif
                else
                  conmap(imap) = con(it,ir,ig)
                endif
                if(imap.gt.2) then
                  conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
                else
                  conmap(1) = conmap(2)
                endif
              enddo  !  it
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=1,nt)
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
              if(monplot) then
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=1,nt)
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=nt+1,2*nt)
              endif
            enddo  !  ir
          enddo  !  ig
        elseif((ng*nr/2).le.(4*nt)) then
          nxm = 2*nt
          nym = ng*nr2/2
          do ig = 1,ng
            do ir = nr2+1,nr,2
              imap = 0
              do it = 1,nt
                imap = imap+2
                if(logplot) then
                  prod = con(it,ir,ig)*con(it,ir+1,ig)
                  if(prod.gt.zero) then
                    conmap(imap) = slog10(prod)/two
                  else
                    conmap(imap) = zero
                  endif
                else
                  conmap(imap) = (con(it,ir,ig)+con(it,ir+1,ig))/two
                endif
                if(imap.gt.2) then
                  conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
                else
                  conmap(1) = conmap(2)
                endif
              enddo  !  it
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=1,nt)
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
              if(monplot) then
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=1,nt)
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=nt+1,2*nt)
              endif
            enddo  !  ir
          enddo  !  ig
        else
          nxm = 2*nt
          nym = ng*nr2/4
          do ig = 1,ng
            do ir = nr2+1,nr,4
              imap = 0
              do it = 1,nt
                imap = imap+2
                if(logplot) then
                  prod = con(it,ir,ig)*con(it,ir+1,ig)
     &                *con(it,ir+2,ig)*con(it,ir+3,ig)
                  if(prod.gt.zero) then
                    conmap(imap) = slog10(prod)/four
                  else
                    conmap(imap) = zero
                  endif
                else
                  conmap(imap) = (con(it,ir,ig)+con(it,ir+1,ig)
     &              +con(it,ir+2,ig)+con(it,ir+3,ig))/four
                endif
                if(imap.gt.2) then
                  conmap(imap-1) = (conmap(imap-2)+conmap(imap))/two
                else
                  conmap(1) = conmap(2)
                endif
              enddo  !  it
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=1,nt)
              irec = irec+1
              write(unit=9,rec=irec) (conmap(it),it=nt+1,2*nt)
              if(monplot) then
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=1,nt)
                jrec = jrec+1
                write(unit=10,rec=jrec) (conmap(it),it=nt+1,2*nt)
              endif
            enddo  !  ir
          enddo  !  ig
        endif

        close(unit=9)
        if(monplot) close(unit=10)
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("contours")')
        write(8,'(3i8)') nxm,nym,1
        close(unit=8)

 90     return
      end



      subroutine plotim(frame,mx,nx,monfile)
c	make greyscale plot of an image frame

        implicit real*8 (a-h,o-z)
        dimension frame(mx,mx)
        character*16 lockline,monfile
        logical monplot

        monplot = (index(monfile,' ').gt.1)
   88   open(unit=8,file='plotlock',access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 88
        endif
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='images',access='direct',recl=2*nx)
        if(monplot) open(unit=10,file=monfile,access='direct',recl=2*nx)

        do iy = 1,nx
          write(unit=9,rec=iy) (frame(ix,iy),ix=1,nx)
          if(monplot) write(unit=10,rec=iy) (frame(ix,iy),ix=1,nx)
        enddo

        close(unit=9)
        if(monplot) close(unit=10)
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("images")')
        write(8,'(3i8)') nx,nx,1
        close(unit=8)

   90   return
      end


      subroutine plotggfs(gg,jg,ju)

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mu = 64)
        parameter (ma = 32)
        parameter (mc = 1024)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (hun = 1.0d+02)

        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        dimension gg(mu,mt,mr,mg)
        dimension frame(mc,mc)
        logical badnum

        nc = mc
        alog2 = log(two)
        call startgr(jg,ig,kr,lr,lt,dr,dcost)
        if(ig.gt.ng) then
          dxm = (nr9*dr9)/(two*nc)
          x0 = (nc-half)*dxm
          zstr = nn9
        else
          dxm = (nr*dr)/nc
          x0 = -half*dxm
          zstr = one
        endif
        rnr = nr*dr0
        summap = zero
        summi = zero
        do ix = 1,nc
          do iz = 1,nc
            x = x0+ix*dxm
            y = zero
            z = (iz-nc/2-half)*dxm
c	stretch z scale
            z = z/zstr
            r = sqrt(x**2+y**2+z**2)
            costh = abs(z/r)
            if(ig.gt.ng) then
              it = nt-nint(costh/dcost-half)
              if((it.ge.1).and.(it.le.nt)) then
                call startgr(jg,kg,kr,lr,lt,drgr,dcosgr)
                ir = r/dr9+one
                ir = min(nr9,max(1,ir))
              else
                igr = slog(r/rnr)/alog2+one
                igr = min(ng,max(1,igr))
                call startgr(igr,kg,kr,lr,lt,drgr,dcosgr)
                ir = r/drgr+one
                it = nt-nint(costh/dcosgr-half)
                if((it.lt.1).or.(it.gt.nt).or.(ir.lt.1).or.(ir.gt.nr))
     &            then
                  frame(ix,iz) = one
                  goto 240
                endif
              endif
            else
              igr = slog(r/rnr)/alog2+one
              igr = min(ng,max(1,igr))
              call startgr(igr,kg,kr,lr,lt,drgr,dcosgr)
              ir = r/drgr+one
              it = nt-nint(costh/dcosgr-half)
              if((it.lt.1).or.(it.gt.nt).or.(ir.lt.1).or.(ir.gt.nr))
     &          then
                frame(ix,iz) = one
                goto 240
              endif
            endif
            if(r.le.rstar) then
              frame(ix,iz) = one
              goto 240
            elseif(r.lt.(hun*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/two
            else
              ostar = (rstar/(two*r))**2
            endif
            if(ju.eq.0) then
              ggsum = zero
              fssum = zero
              do iu = 1,nu
                ggsum = ggsum+opa(iu,nt,nr,1)*gg(iu,it,ir,kg)*dud(iu)
                fssum = fssum+opa(iu,nt,nr,1)*bus(iu)*dud(iu)
              enddo
              if(fssum.gt.zero) then
                ggfs = ggsum/(fssum*ostar)
              else
                ggfs = one
              endif
            elseif(ju.lt.0) then
              ia = abs(ju)
              ggsum = zero
              fssum = zero
              do iu = 1,nu
                ggsum = ggsum+phab(iu,ia)*gg(iu,it,ir,kg)*dud(iu)
                fssum = fssum+phab(iu,ia)*bus(iu)*dud(iu)
              enddo
              if(fssum.gt.zero) then
                ggfs = ggsum/(fssum*ostar)
              else
                ggfs = one
              endif
            else
              ggfs = gg(ju,it,ir,kg)/(bus(ju)*ostar)
            endif
            if(badnum(ggfs)) then
              print '(" in ggfs:",5i6,1p,2e10.2)',
     &          ix,iz,it,ir,kg,gg(ju,it,ir,kg),bus(iu)
              ggfs = zero
            endif
            summap = summap+ggfs
            summi = summi+one
            frame(ix,iz) = sqrt(ggfs)
  240     enddo
        enddo

        print '(" plotting sqrt(gg/stellar) with z stretch =",f4.1)',
     &    zstr
        print '(" mean gg/fs =",1p,e10.3)', summap/summi
        call plotim(frame,mc,nc,'')

        return
      end


      subroutine plotxyim(pl,jg)

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (mc = 1024)
        parameter (zero = 0.0d+00)
        parameter (half = 0.5d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)

        common /mpis/mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        dimension pl(mt,mr,mg)
        dimension frame(mc,mc)

        nc = mc
        alog2 = log(two)
        call startgr(jg,ig,kr,lr,lt,dr,dcost)
        if(ig.gt.ng) then
          rnr = kr*dr
          dxm = rnr/(2*nc)
          x0 = (nc-half)*dxm
        else
          dxm = dr
          rnr = nr*dr0
          x0 = -half*dxm
        endif
        summap = zero
        summi = zero
        do ix = 1,nc
          do iz = 1,nc
            x = x0+ix*dxm
            y = zero
            z = (iz-nc/2-half)*dxm
            r = sqrt(x**2+y**2+z**2)
            costh = abs(z/r)
            if(ig.gt.ng) then
              it = nt-nint(costh/dcost-half)
              if((it.ge.1).and.(it.le.nt)) then
                igr = ig
                drgr = dr9
                ir = r/dr9+one
                ir = min(nr9,max(1,ir))
              else
                igr = slog(r/rnr)/alog2+one
                igr = min(ng,max(1,igr))
                call startgr(igr,kg,kr,lr,lt,dr,dcosgr)
                drgr = dr
                ir = r/drgr+one
                ir = min(nr,max(1,ir))
                it = nt-nint(costh/dcosgr-half)
                if((it.lt.1).or.(it.gt.nt)) then
                  frame(ix,iy) = zero
                  goto 240
                endif
              endif
            else
              igr = slog(r/rnr)/alog2+one
              igr = min(ng,max(1,igr))
              call startgr(igr,kg,kr,lr,lt,dr,dcosgr)
              drgr = dr
              ir = r/drgr+one
              ir = min(nr,max(1,ir))
              it = nt-nint(costh/dcosgr-half)
              if((it.lt.1).or.(it.gt.nt)) then
                frame(ix,iy) = zero
                goto 240
              endif
            endif
            if(r.le.rstar) then
              frame(ix,iz) = zero
              goto 240
            elseif(r.lt.(100.*rstar)) then
              ostar = (one-sqrt(one-(rstar/r)**2))/2.
            else
              ostar = 0.25*(rstar/r)**2
            endif
c            frame(ix,iz) = pl(it,ir,kg)/ostar
            frame(ix,iz) = pl(it,ir,kg)*(r/rstar)
            summap = summap+frame(ix,iz)
            summi = summi+one
            frame(ix,iz) = max(frame(ix,iz),zero)
  240     enddo
        enddo

        call plotim(frame,mc,nc,'')

        return
      end


      subroutine plotspec(spec,mw,nw,nspec,irot,dw,monfile)
c	plot spectrum

        implicit real*8 (a-h,o-z)
        dimension spec(mw,nspec)
        character*16 lockline,monfile
        logical monplot

        monplot = (index(monfile,' ').gt.1)
   88   open(unit=8,file='plotlock',access='sequential',err=90)
        read(8,'(a16)') lockline
        close(unit=8)
        if(index(lockline,'reading').gt.0) then
          call sleep(1)
          goto 88
        endif
        open(unit=8,file='plotlock',access='sequential',err=90)
        write(8,'("writing")')
        close(unit=8)
        open(unit=9,file='specs',access='direct',recl=2*nw)
        if(monplot) open(unit=10,file=monfile,access='direct',recl=2*nw)

        do ispec = 1,nspec
          write(unit=9,rec=ispec) (spec(iw,ispec),iw=1,nw)
          jspec = (irot-1)*nspec+ispec
          if(monplot) write(unit=10,rec=jspec) (spec(iw,ispec),iw=1,nw)
        enddo

        close(unit=9)
        if(monplot) close(unit=10)
        open(unit=8,file='plotlock',access='sequential')
        write(8,'("specs")')
        write(8,'(2i8,1p,e12.3)') nw,nspec,dw
        close(unit=8)

   90   return
      end



      function bw(w,t)
c	black body function

        implicit real*8 (a-h,o-z)
        parameter (hc2 = 5.955d-06)
        parameter (hck = 1.4388d+00)
        parameter (big = 5.0d+02)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)

        arg = hck*w/t
        if(arg.le.zero) then
          bw = zero
        elseif(arg.lt.big) then
          bw = two*hc2*w**3/(sexp(arg)-one)
        else
          bw = zero
        endif

        return
      end


      function texp(arg)
c	exp(arg) using look-up table

        implicit real*8 (a-h,o-z)
        common /exptab/ expi(10239),too
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (hun = 1.0d+02)
        parameter (big = 1.0d+04)

        harg = abs(hun*arg)
        if(harg.lt.big) then
          iarg = int(harg)
          darg = (harg-iarg)/hun
          texp = (one-darg+darg**2/too)*expi(iarg+1)
          if(arg.gt.zero) texp = one/texp
        elseif(arg.lt.zero) then
          texp = zero
        else
          texp = one/expi(10000)
        endif

        return
 
      entry expinit

        do iarg = 1,10000
          expi(iarg) = sexp(-(iarg-1)/hun)
        enddo
        darg = one/hun
        too = darg**2/(exp(-darg)-one+darg)
        expinit = one

        return
      end


      function sexp(arg)
c	safe and efficient exp function

        implicit real*8 (a-h,o-z)
        parameter (big = 4.6d+02)
        parameter (small = 1.0d-03)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        save nan
        if((nan.lt.0).or.(nan.gt.8)) nan = 0

        if(arg.lt.(-big)) then
          sexp = zero
        elseif(arg.gt.big) then
          sexp = 1.0d+200
        elseif(abs(arg).lt.small) then
          sexp = 1.0+arg+arg**2/two
        elseif((arg.lt.zero).or.(arg.gt.zero)) then
          sexp = exp(arg)
        else  !  NaN
          print '(e10.2," passed to sexp",i6)', arg,nan
          nan = nan+1
          if(nan.gt.8) stop
          sexp = one
        endif

        return
      end


      function slog(arg)
c	safe log function

        implicit real*8 (a-h,o-z)
        parameter (tiny = 1.0d-200)
        parameter (big = 1.0d+200)
        parameter (small = 1.0d-03)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        save nan
        if((nan.lt.0).or.(nan.gt.8)) nan = 0

        if(arg.le.zero) then
          slog = -500.
        elseif(arg.lt.tiny) then
          slog = -460.
        elseif(arg.gt.big) then
          slog = 460.
        elseif(abs(arg-one).lt.small) then
          slog = arg-one
          slog = slog-slog**2/two
        elseif(arg.gt.tiny) then
          slog = log(arg)
        else  !  NaN
          if(nan.le.8) print '(e10.2," passed to slog",i6)', arg
          nan = nan+1
          slog = zero
        endif

        return
      end


      function slog10(arg)
c	safe log_10 function

        implicit real*8 (a-h,o-z)
        parameter (tiny = 1.0d-200)
        parameter (big = 1.0d+200)
        parameter (small = 1.0d-03)
        parameter (aloge = 0.4342945d+00)
        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        save nan
        if((nan.lt.0).or.(nan.gt.8)) nan = 0

        if(arg.le.zero) then
          slog10 = -250.
        elseif(arg.lt.tiny) then
          slog10 = -200.
        elseif(arg.gt.big) then
          slog10 = 200.
        elseif(abs(arg-one).lt.small) then
          slog = arg-one
          slog10 = aloge*(slog-slog**2/two)
        elseif(arg.gt.tiny) then
          slog10 = log10(arg)
        else  !  NaN
          if(nan.le.8) print '(e10.2," passed to slog10",i6)', arg
          nan = nan+1
          slog10 = zero
        endif

        return
      end


      subroutine waitasec(wait,iwait)

        implicit real*8 (a-h,o-z)
        common /mpis/mysize,myrank
        logical wait
        character*16 lockline

        if(wait) then
          call sleep(iwait-1)
  100     call sleep(1)
          open(unit=8,file='plotlock',access='sequential',err=190)
          read(8,'(a16)') lockline
          close(unit=8)
          if(index(lockline,'empty').eq.0) goto 100
        else
          call sleep(iwait)
        endif

  190   return
      end


      subroutine savemodel

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)


        open(unit=9,file='diskmodel',access='direct',recl=2*nt)
        irec = 0
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            irec = irec+1
            write(unit=9,rec=irec) (dg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            irec = irec+1
            write(unit=9,rec=irec) (tg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            irec = irec+1
            write(unit=9,rec=irec) (opl(ml,it,ir,ig),it=1,nt)
          enddo
        enddo
        close(unit=9)
        irec = 0
        open(unit=9,file='dustmodel',access='direct',recl=2*nu)
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,nt
              irec = irec+1
              write(unit=9,rec=irec) (opa(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,nt
              irec = irec+1
              write(unit=9,rec=irec) (swd(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        close(unit=9)

        return
      end


      subroutine readmodel

        implicit real*8 (a-h,o-z)
        parameter (mg = 20)
        parameter (mr = 128)
        parameter (mt = 192)
        parameter (ma = 32)
        parameter (ml = 2)
        parameter (mu = 64)
        parameter (mw = 1024)
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /dust/ ad(ma),popd(ma),wd(mu),ud(mu),dud(mu),bus(mu),
     &          phid(mu,ma),phab(mu,ma),phsc(mu,ma),dd(mt,mr,mg),
     &          opa(mu,mt,mr,mg),sca(mu,mt,mr,mg),swd(mu,mt,mr,mg),
     &          emis(mu),rebus(mu),condens(ma,mt,mr,mg),td(mt,mr,mg)
        common /gas/ dg(mt,mr,mg),tg(mt,mr,mg),gctg(mt,mr,mg),fstar(ml),
     &          wl(ml),opl(ml,mt,mr,mg)
        logical badnum


        open(unit=9,file='diskmodel',access='direct',recl=2*nt,err=98)
        irec = 0
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (dg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (tg(it,ir,ig),it=1,nt)
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            irec = irec+1
            read(unit=9,rec=irec,err=98) (opl(ml,it,ir,ig),it=1,nt)
          enddo
        enddo
        close(unit=9)

        irec = 0
        open(unit=9,file='dustmodel',access='direct',recl=2*nu,err=99)
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,nt
              irec = irec+1
              read(unit=9,rec=irec,err=99) (opa(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,nt
              irec = irec+1
              read(unit=9,rec=irec,err=99) (swd(iu,it,ir,ig),iu=1,nu)
            enddo
          enddo
        enddo
        close(unit=9)

        do igr = 0,ng
          call startgr(igr,ig,kr,lr,lt,dr,dcosth)
          do ir = 1,lr
            do it = 1,nt
              if(badnum(dg(it,ir,ig)))
     7          print '("dg(",3i4,") =",e10.2)', it,ir,ig,dg(it,ir,ig)
              if(badnum(tg(it,ir,ig)))
     7          print '("tg(",3i4,") =",e10.2)', it,ir,ig,tg(it,ir,ig)
              if(badnum(opl(ml,it,ir,ig))) print
     7          '("opl(",3i4,") =",e10.2)', it,ir,ig,opl(ml,it,ir,ig)
              do iu = 1,nu
                if(badnum(opa(iu,it,ir,ig)))
     7            print '("opa(",4i4,") =",e10.2)',
     &            iu,it,ir,ig,opa(iu,it,ir,ig)
                if(badnum(swd(iu,it,ir,ig)))
     7            print '("swd(",4i4,") =",e10.2)',
     &            iu,it,ir,ig,swd(iu,it,ir,ig)
              enddo
            enddo
          enddo
        enddo

        return

   98   print '(" Error reading diskmodel")'
        stop

   99   print '(" Error reading dustmodel")'
        stop

      end



      function badnum(val)
c       test val for NaN or Inf

        implicit real*8 (a-h,o-z)
        logical badnum

        if(val.eq.(0.0d+00)) then
          badnum = .false.
        else
          test = val/val
          badnum = (test.ne.test)
        endif

        return
      end



      subroutine readparms

        implicit real*8 (a-h,o-z)
        character*80 line,parn,parv
        common /mpis/ mysize,myrank
        common /nn/ nx,ny,nz,ng,nw,nr,nt,na,nu,nl,nrot,lw,ll,nit
        common /nines/ nx9,nz9,nr9,ng9,nn9,mm9,dx9,dz9,dr9
        common /consts/ we,be,ae,hc,sv,wec,acoeff,gnucl,popiso,moltype
        common /parms/ dx0,dz0,dr0,dv,dvturb0,dvturb1,vterm,rstar,tstar
        common /heat/ gmmdot,acom,pxr,axr,epspe,agas,sigma,
     &          colin,colup,xco,xwat
        common /angls/ dcost,dcos9,rot,irot
        common /parmz/ rgrid1,r0,sigma0,pcol,qcol,avisc,xmstar,xlstar,
     &          tuv,xuv,tsubl,gsubl,enhanc,tdiss,tice,xco0,xwat0,
     &          coswind,xc,xsi,ac,asi,kit,lit,mit,nnt,lx

        parameter (zero = 0.0d+00)
        parameter (one = 1.0d+00)
        parameter (two = 2.0d+00)
        parameter (half = 0.5d+00)
        parameter (quart = 0.25d+00)
        parameter (pi = 3.14159265d+00)
        parameter (au = 1.496d+13)
        parameter (rsun = 6.96d+10)
        parameter (xlsun = 3.90d+33)
        parameter (xmsun = 1.99d+33)
c       mass / gas particle
        parameter (xmg = 4.0d-24)
        parameter (vkms = 1.0d+05)

        open(unit=12,file='disco.par')

  100   read(unit=12,fmt='(a72)',err=90,end=90) line
        if(index(line,'c').eq.1) then
          continue
        elseif(index(line,'=').gt.0) then
          npar = index(line,'=')-1
          parn = line(1:npar)
          npar = npar+2
          parv = line(npar:72)
          if(index(parn,'rgrid1').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) rgrid1
            rgrid1 = rgrid1*au
          elseif(index(parn,'r0').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) r0
            r0 = r0*au
          elseif(index(parn,'sigma0').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) sigma0
            sigma0 = sigma0/xmg
          elseif(index(parn,'pcol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) pcol
          elseif(index(parn,'qcol').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) qcol
          elseif(index(parn,'avisc').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) avisc
            avisc = avisc*xmg/3.0d+00
          elseif(index(parn,'rstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) rstar
            rstar = rstar*rsun
          elseif(index(parn,'mstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xmstar
          elseif(index(parn,'lstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xlstar
            tstar = 5.8d+03*(xlstar/(rstar/rsun)**2)**quart
          elseif(index(parn,'tstar').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tstar
            xlstar = (rstar/rsun)**2*(tstar/5.8d+03)**4
          elseif(index(parn,'tuv').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tuv
          elseif(index(parn,'xuv').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xuv
          elseif(index(parn,'tsubl').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tsubl
          elseif(index(parn,'gsubl').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) gsubl
          elseif(index(parn,'enhanc').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) enhanc
          elseif(index(parn,'tdiss').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tdiss
          elseif(index(parn,'tice').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) tice
          elseif(index(parn,'mdot').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xmdot
          elseif(index(parn,'acom').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) acom
          elseif(index(parn,'pxr').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) pxr
          elseif(index(parn,'axr').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) axr
          elseif(index(parn,'epspe').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) epspe
          elseif(index(parn,'xco').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xco0
          elseif(index(parn,'xwat').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xwat0
          elseif(index(parn,'coswind').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) coswind
          elseif(index(parn,'dvturb0').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) dvturb0
            dvturb0 = dvturb0*vkms
          elseif(index(parn,'dvturb1').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) dvturb1
            dvturb1 = dvturb1*vkms
          elseif(index(parn,'dv').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) dv
            dv = dv*vkms
          elseif(index(parn,'xc').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xc
          elseif(index(parn,'xsi').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) xsi
          elseif(index(parn,'ac').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) ac
          elseif(index(parn,'asi').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) asi
          elseif(index(parn,'kit').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) kit
          elseif(index(parn,'lit').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) lit
          elseif(index(parn,'mit').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) mit
          elseif(index(parn,'ng').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) ng
          elseif(index(parn,'nn9').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nn9
          elseif(index(parn,'nr').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nr
          elseif(index(parn,'nnt').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nnt
          elseif(index(parn,'nt').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nt
          elseif(index(parn,'nx').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) nx
          elseif(index(parn,'lx').gt.1) then
            read(unit=parv,fmt=*,err=90,end=90) lx
          else
            print '("I don''t understand",a60)', line(1:60)
          endif
        endif
        goto 100

   90   close(unit=12)
        return
      end
